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This work considers a theoretical model for quasi-spherical subsonic accretion onto slowly rotating mag- 
, netized neutron stars. In this regime the accreting matter settles down subsonically onto the rotating 

magnetosphere, forming an extended quasi-static shell. The shell mediates the angular momentum trans- 
\ for to/from the rotating neutron star magnetosphere by large-scale convective motions, which for observed 

■ pulsars lead to an almost so-angular-momentum rotation law with w ~ inside the shell. The ac- 

, cretion rate through the shell is determined by the ability of the plasma to enter the magnetosphere due 

^ ■ to Ray leigh- Taylor instabilities while taking cooling into account. The settling regime of accretion is 

possible for moderate accretion rates M < M^, ~ 4 x 10"'^^ g/s. At higher accretion rates a free-fall gap 
CO ' above the neutron star magnetosphere appears due to rapid Compton cooling, and accretion becomes 

highly non-stationary. From observations of spin- up/spin-down rates of quasi-spherically wind accreting 
j-j"! '■ equilibrium X-ray pulsars with known orbital periods (like e.g. GX 301-2 and Vela X-1), it is possible to 

determine the main dimensionless parameters of the model, as well as to estimate the magnetic field on 
the surface of the neutron star. For equilibrium pulsars with independent measurements of the magnetic 
I field, the model also allows us to estimate the velocity of the stellar wind from the companion without the 

use of complicated spectroscopic measurements. For non-equilibrium pulsars, it can be shown that there 



Q , exists a maximum possible value of the spin-down rate of the accreting neutron star. From observations of 

the spin-down rate and the X-ray luminosity in such pulsars (e. g. GX 1-1-4, SXP 1062 and 4U 2206-f 54) 
^ I we are able to estimate a lower limit on the neutron star magnetic field, which in all exemplified cases 

turns out to be close to the standard one and in agreement with cyclotron line measurements. The model 
further explains both the spin- up/spin-down of the pulsar frequency on large time-scales and the irreg- 
. ular short-term frequency fiuctuations, which may correlate or anti-correlate with the X-ray luminosity 

I fiuctuations, seen in different systems. 

o 

O ; 1 Introduction 

' X-ray pulsars are highly magnetized rotating neutron stars in close binary systems, accreting matter 

^■f-^ . from a companion star. The companion may be a low-mass star overfilling its Roche lobe, in which case 

an accretion disc is formed. In the case of a high-mass companion, the neutron star may also accrete 
from the strong stellar wind, and depending on the conditions a disc may be formed or accretion may 
take place quasi-spherically. The strong magnetic field (of the order of 10^^ — 10^^ G) of the neutron 
. star disrupts the accretion fiow at some distance from the neutron star surface and forces the accreted 

matter to funnel down on the polar caps of the neutron star creating hot spots that, if misaligned with 
the rotational axis, make the neutron star pulsate in X-rays. Most accreting pulsars show stochastic 
variations in their spin frequencies as well as in their luminosities. Many sources also exhibit long-term 
trends in their spin-behaviour with the period more or less steadily increasing or decreasing, and in some 
sources spin-reversals have been observed. (For a thorough review, see e.g. [1] and references therein.) 

The best-studied case of accretion is that of thin disc accretion [2]. Here the spin- up/spin-down 
mechanisms are rather well understood. For disc accretion the spin-up torque is determined by the specific 
angular momentum at the inner edge of the disc and can be written in the form K^u ~ M\/GMRa [3j . 
For a pulsar the inner radius of the accretion disc is determined by the Alfven radius Ra ^ M~^/^, 
so Ksu ^ M^f , i.e. for disc accretion the spin- up torque is weakly (almost linearly) dependent on 
the accretion rate (X-ray luminosity). In contrast, the spin-down torque for disc accretion in the first 
approximation is independent of M: Ksd ^ /Rci where Rc = (G'M/(w*)^)^^'^ is the corotation radius, 
uj* is the neutron star angular frequency and fi is the neutron star's dipole magnetic moment. In fact, 
accretion torques in disc accretion are determined by complicated disc-magnetospheric interactions, see, 
6-g-, and the discussion in [6;, and correspondingly can have a more complicated dependence on 

the mass accretion rate and other parameters. 

Measurements of spin- up/spin-down in X-ray pulsars can be used to evaluate a very important param- 
eter of the neutron star - its magnetic field. The period of the pulsar is usually close to the equilibrium 
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value Peq, which is determined by the total zero torque applied to the neutron star, K ~ Kg^ + K^d ~ 0. 
So assuming the observed value w* — 2t: / P^q^ the magnetic field of the neutron star in disc-accreting 
X-ray pulsars can be estimated if M is known. 

In the case of quasi-spherical accretion, which may take place in systems where the optical star 
underfills its Roche lobe and no accretion disc is formed, the situation is more complicated. Clearly, 
the amount and sign of the angular momentum supplied to the neutron star from the captured stellar 
wind are important for spin-up or spin-down. To within a numerical factor of the order of 1 (which can 
be positive or negative, see numerical simulations by [2], [8], [9], etc.), the torque applied to the neutron 
star in this case should be proportional to MujbR%, where ujb = 2'k/Pb is the binary orbital angular 
frequency, Rb — 2GM/{V^ + vl^^^Y is the gravitational capture (Bondi) radius, is the steUar wind 
velocity at the neutron star orbital distance, and Vorb is the neutron star orbital velocity. In real high-mass 
X-ray binaries the orbital eccentricity is non-zero, the stellar wind is variable and can be inhomogeneous, 
etc., so can be a complicated function of time. The spin-down torque is even more uncertain, since 
it is impossible to write down a simple equation like —^^/R^ any more {Rc has no meaning for quasi- 
spherical accretion; for slowly rotating pulsars it is much larger than the Alfven radius where the angular 
momentum transfer from the accreting matter to the magnetosphere actually occurs). For example, 
using the expression —ji^/R^ for the braking torque results in a very high (> 10^'' G) magnetic field for 
long-period X-ray pulsars. We think this is a result of underestimating the braking torque. 

The matter captured from the stellar wind can accrete onto the neutron star in different ways. Indeed, 
if the X-ray flux from the accreting neutron star is sufficiently high, the shocked matter rapidly cools down 
due to Compton processes and falls freely toward the magnetosphere. The velocity of motion rapidly 
becomes supersonic, so a shock is formed above the magnetosphere. This regime was considered, e.g., by 
[To] . Depending on the sign of the specific angular momentum of falling matter (prograde or retrograde), 
the neutron star can spin-up or spin-down. However, if the X-ray flux at the Bondi radius is below 
some value, the shocked matter remains hot, the radial velocity of the plasma is subsonic, and the source 
may enter the settling accretion regime. A hot quasi-static shell forms around the magnetosphere [TT] . 
Due to additional energy release (especially near the base of the shell), the temperature gradient across 
the shell becomes superadiabatic, and large-scale convective motions inevitably appear. The convection 
initiates turbulence, and the motion of a fluid element in the shell becomes quite complicated. If the 
magnetosphere allows plasma entry via instabilities (and subsequent accretion onto the neutron star), 
the actual accretion rate through such a shell is controlled by the magnetosphere (for example, a shell 
can exist, but accretion through it can be weak or even absent altogether). Therefore, on top of the 
convective motions, the matter acquires a low, on average radial, velocity toward the magnetosphere, 
and thus subsonic settling is possible. This type of accretion can work only for small X-ray luminosities. 
La; < 4 X 10'^^ erg/s (see below), and is totally different from that considered in the numerical simulations 
cited above. If a shell is present, its interaction with the rotating magnetosphere can lead to spin-up or 
spin-down of the neutron star, depending on the sign of the difference of the angular velocity between the 
accreting matter and the magnetospheric boundary. Thus, in the settling accretion regime, both spin-up 
or spin-down of the neutron star is possible, even if the sign of the specific angular momentum of the 
captured matter is always prograde. The shell here mediates the angular momentum transfer to or from 
the rotating neutron star. 

There are several models in the literature (see especially and H^), from which the expression for 
the spin-down torque for quasi-spherically accreting neutron stars in the form K^d ~ —MR^lo* can be 
derived. Using the standard expression for the Alfven radius, Ra ^ Af ^^/^/i'*/'' this torque is proportional 
to Ksd ~ —[^I'^M'^I'^. In our model, the matter in the shell settles subsonically as the region close to 
the magnetospheric surface cools down, and the Alfven radius has a different dependence on the mass 
accretion rate and the magnetic field, Ra ^ A/^^/^^/^^/^^ (see below). 

One can show that there are two different mechanisms through which angular momentum can be 
transferred through a quasi-spherical shell. In the first case (we call this case moderate coupling), angular 
momentum is transferred by convective motions in the shell. The breaking torque in the regime of settling 
accretion with convective removal of angular momentum depends on the accretion rate as K^d ~ —M'^l^'^ 
(see Section HI). The velocity of the convective motions in this regime is close to the sound speed. It is also 
possible to have a settling regime where the angular momentum is removed by shear turbulence in the 
shell (the weak coupling regime). In this regime the characteristic velocities of the shear fiow close to the 
magnetosphere is of the order of the linear rotational velocity. In this case K^d ~ lR?c ^ M^i^*^/ {GM\ 
i.e. in the weak coupling regime the torque does not depend on the accretion rate at all. 
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Figure 1: A schematic picture of quasi-spherical accretion from the stellar wind of an optical companion 
star (left) onto a neutron star (right) in a binary system. In the regime of subsonic accretion, a quasi- 
spherical shell (shaded area) of radius Ra is formed between the bow shock (parabolic curve) and the 
rotating magnetosphere. In this shell, large-scale convective motions are formed that may act to remove 
angular momentum from the magnetosphere. The outer radius of the shell is determined by the gravi- 
tational capture or Bondi radius Rb- The characteristic velocity of the wind is ~ 300 — 1000 km/s. 
The contour arrow shows the direction of the orbital velocity of the neutron star Vorb- 
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To stress the difference between the two possible regimes of subsonic accretion (with moderate 
and weak couphng), we can rewrite the expression for the breaking torque with convection (moder- 
ate couphng) using the corotational radius and the Alfven radius in the form Kgd ~ —fJ-^ / VRc^a ~ 
-{n"^ / Rl){Rc/ RAf'"^ (see further details in Section 3). Since the factor {Rc/RaY^'^ ^ {uJKiRA)/^^*) can 
be of the order of 10 or more in real systems, using a braking torque in the form of may lead to a 

strong overestimate of the magnetic field of the neutron star. 

The dependence of the braking torque on the accretion rate in the case of quasi-spherical settling 
accretion suggests that variations of the mass accretion rate (and X-ray luminosity) must lead to a 
transition from spin-up (at high accretion rates) to spin-down (at small accretion rates) at some critical 
value of M (or Ra), that differs from source to source. This phenomenon (known as torque reversal) 
is actually observed in wind-fed pulsars like Vela X-1, GX 301-2 and GX l-f-4, which we shall consider 
below in more detail. 

The structure of this paper is as follows. In Section 2, we present an outline of the theory for quasi- 
spherical accretion onto a neutron star magnetosphere. We show that it is possible to construct a hot 
envelope around the neutron star through which subsonic accretion can take place and act to either spin 
up or spin down the neutron star. In Section 3, we discuss the structure of the interchange instability 
region which determines whether the plasma can enter the magnetosphere of the rotating neutron star. In 
Section 4 we consider how the spin- up/spin-down torques vary with a changing accretion rate. In Section 
5, we show how to determine the parameters of quasi-spherical accretion from observational data. In 
Section 6, we apply our methods to the specific pulsars GX 301-2, Vela X-1, GX 1+4, SXP 1062 and 
4U 2206+54. In Section 7 we discuss our results and, finally, in Section 8 we present our conclusions. A 
detailed gas-dynamic treatment of the problem is presented in five appendices, which are very important 
to understand the physical processes involved. 

This work follows to a large extent the earlier published paper of |14) . However, here are included 
several additions, clarifying and refining the physical model (especially in Sections 1-4 and in Conclusions. 

2 Quasi-spherical accretion 

2.1 The structure of a subsonic shell around a neutron star magnetosphere 

We shall here consider the torques applied to a neutron star in the case of quasi-spherical accretion from 
a stellar wind. Wind matter is gravitationally captured by the moving neutron star and a bow-shock 
is formed at a characteristic distance R Rb, where Rb is the Bondi radius. Angular momentum can 
be removed from the neutron star magnetosphere in two ways — either with matter expelled from the 
magnetospheric boundary without accretion (the propeller regime, |15j). or via large-scale convective 
motions in a subsonic quasi-static shell around the magnetosphere, in which case the accretion rate onto 
the neutron star is determined by the ability of the plasma to enter the magnetosphere, in the regime of 
subsonic accretion. 

In such a quasi-static shell, the temperature will be high (of the order of the virial temperature, see 
|11|). and the important point is whether hot matter from the shell can in fact enter the magnetosphere. 
Two-dimensional calculations by [TB] have shown that hot monoatomic ideal plasma is stable relative 
to the Rayleigh- Taylor instability at the magnetospheric boundary, and plasma cooling is thus needed 
for accretion to begin. However, a closer inspection of the 3-dimensional calculations by [17] reveals 
that the hot plasma is only marginally stable at the magnetospheric equator (to within 5% accuracy 
of their calculations). Compton cooling and the possible presence of dissipative phenomena (magnetic 
reconnection etc.) facilitates the plasma entering the magnetosphere. We will show below that spin- 
down of the neutron star is possible in the case of accretion of matter from a hot envelope in the subsonic 
settling regime. 

To a zeroth approximation, we can neglect both rotation and radial motion (accretion) of matter in 
the shell and consider only its equilibrium hydrostatic structure. The radial velocity of matter falling 
through the shell un is lower than the sound velocity Cg- Under these assumptions, the characteristic 
cooling/heating time-scale is much larger than the free-fall time-scale. 

In the general case where both gas pressure and anisotropic turbulent motions are present, Pascal's 
law is violated. Then the hydrostatic equilibrium equation can be derived from the equation of motion 
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(|A.16p with stress tensor components (IA.19P - (|A.2ip and zero viscosity (see Appendix A for more detail) : 



IdP, 1 diPl^R"^) 2P\ GM ^ 

1 — =0 fl) 

p dR pR"^ dR pR R^ ^ ' 

Here Pg ~ pc^sll the gas pressure, and P* stands for the pressure due to turbulent motions: 

P|f = p < u\ pmjjc^ = 7Pgmjj (2) 

Pi^p< u\ >= pmlcl - ^Pgml (3) 

(< >=< u| > +2 < > is the turbulent velocity dispersion, and are turbulent Mach 
numbers squared in the radial and tangential directions, respectively; for example, in the case of isotropic 
turbulence to| = nr^ = (l/3)m( where rrit is the turbulent Mach number). The total pressure is the sum 
of the gas and turbulence terms: Pg + Pt ~ Pg{l + ^nif). 

The turbulent Mach number in the shell may in general depend on the radius. In our case, however, 
we will consider it constant. Furthermore, in real pulsars, turbulent heating (important from a dynamic 
point of view, see Appendix E) will change the estimated parameters by less than a factor of 2 (see 
formulas in Section 6). 

We shall consider, to a first approximation, that the entropy S is constant throughout the shell. For 
an ideal gas with adiabatic index 7 and equation of state P = Ke^/^"^ p^ , the density can be expressed as 
a function of temperature: p ^ T^/i"/^^), Integrating the hydrostatic equilibrium equation ([T}, we find: 

UT /-f-l\GMf 1 \ 7-IGM,, 



Mm V 7 y \^l + 7w|-2(7-l)(m|-m2_)y 7 R 

(In this solution we have neglected the integration constant, which is not important deep inside the shell. 
It is important in the outer part of the shell, but since the outer region close to the bow shock at ^ i?^ is 
not spherically symmetric, its structure can only be found numerically). We note that taking turbulence 
into account somewhat decreases the temperature within the shell. Most important, however, is that the 
anisotropy of turbulent motions, caused by convection, in the stationary case changes the distribution of 
angular velocity in the shell. Below we will show that in the case of isotropic turbulence, the angular 
velocity distribution within the shell is close to quasi-Keplerian: a;(i?) ^ R^^^^. In the case of strongly 
anisotropic turbulence caused by convection, m| ^ , the distribution of momentum in the shell may 

become almost iso-angular: a;(i?) ^ R^^. Below we shall see that an analysis of several real X-ray pulsars 
favors an iso-angular momentum rotation distribution. 

Now, let us write down how the density varies inside the quasi-static shell for R <C Rb- For a fully 
ionized gas with 7 = 5/3 we find: 

piR) = piR^) (^^y (5) 
Ka 



P(i?) = P{Ra) [-f j . (6) 

The above equations describe the structure of an ideal static adiabatic shell above the magnetosphere. 
Of course, at i? ^ Rb the problem is essentially non-spherically symmetric and numerical simulations 
are required. 

Corrections to the adiabatic temperature gradient due to convective energy transport through the 
shell are calculated in Appendix D. 

2.2 The Alfven surface 

At the magnetospheric boundary (the Alfven surface), the total pressure (including isotropic gas pressure 
and the possibly anisotropic turbulent pressure) is balanced by the magnetic pressure B'^/{8tt) 

Pg + Pt^ P,{Ra){1 + Wt) = ■ (7) 
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The magnetic field at the Alfven radius is determined by the dipole moment and magnetic field of the 
neutron star and by electric currents flowing on the Alfvenic surface (in the magnetopause): 



P,,«,, . ^1^^ |iy . (8) 

(1 + ^mf) Stt \RaJ Mm 

where the dimensionless coefficient K2 takes into account the contribution from these currents and the 
factor 1/(1 + 7m() is due to the turbulent pressure term. For example, in the model by Arons and Lea 
[T7] (their Eq. 31), K2 = (2.75)^ ~ 7.56. At the magnetospheric cusp (where the magnetic force line is 
branched), the radius of the Alfven surface is about 0.51 times that of the equatorial radius [T7]. Below 
we shall assume that Ra is the equatorial radius of the magnetosphere, unless stated otherwise. 

The plasma is able to enter the magnetosphere mainly due to the interchange instability. In the 
stationary regime, let us introduce the accretion rate M onto the neutron star surface. From the continuity 
equation in the shell we find 

M 

47rufl(i?A)i?l 

Clearly, the velocity of absorption of matter by the magnetosphere is smaller than the free-fall velocity, so 
we introduce a dimensionless factor f{u) = urJ yj2GM/R < 1. Then the density at the magnetospheric 
boundary is 

piRA) - , ^, . ^, . (10) 



AiTf{u)^2GM/RAR 



For example, in the model calculations by [T^, f{u) ~ 0.1; in our case, at high X-ray luminosities, the 
value of f{u) may attain « 0.5. If we imagine that the shell is impenetrable and that there is no accretion 
through it, M — )• 0. In this case ur — )■ 0, f{u) — 0, while the density in the shell remains finite. In some 
sense, the matter is leaking from the magnetosphere down onto the neutron star, and the leakage may 
be either very small (M — 0) or have a finite non-zero value (M 7^ 0). 

Plugging p{R) into ([8|) and using (jH) and the definition of the dipole magnetic moment 

M = 2^o^o 

(where Rq is the neutron star radius), we find an expression for the Alfven radius in the case of quasi- 
spherical accretion: 



Ra = 



47 f{u)K2 



(11) 



(7-1) V'(7, "^t)(l + im'^t) mV2GM_ 
It should be stressed that in the presence of a hot shell the Alfven radius is determined by the static 
gas pressure (with a possible addition of turbulent motions) at the magnetospheric boundary, which is 
non-zero even for a zero- mass accretion rate through the shell. The dependence of f{u) on the accretion 
rate M in the case of a settling shell taking cooling into account will be derived below (see ([55)) below) . 
In the supersonic (Bondi) regime we obviously have f{u) = 1. We note that accretion with subsonic 
velocity can take place even in the Bondi regime, but with significantly lower accretion rate (as compared 
to the maximum). In the Bondi regime (i.e. in the adiabatic regime without gas heating and/or cooling), 
the choice of solution depends on the boundary conditions. 

2.3 The mean velocity of matter entering through the magnetospheric bound- 
ary 

As mentioned above, the plasma enters the magnetosphere of the slowly rotating neutron star due to 
the interchange instability. The boundary between the plasma and the magnetosphere is stable at high 
temperatures T > Tcr, but becomes unstable at T < Tcr, and remains in a neutral equilibrium at T = Tcr 
[16] . The critical temperature is: 

■7?T _ 1 COSXf^raGM 

"-'^"^ - 2{l + ^ml) kRa Ra ' ^ ' 

Here k is the local curvature of the magnetosphere, x is the angle the outer normal makes with the 
radius- vector at a given point, and the contribution of turbulent pulsations in the plasma to the total 
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pressure is taken into account by the factor (1 + jm^). The effective gravity acceleration can be written 
as 

GM ( T 



5e// = ^cosx(^l-— j . (13) 

The temperature in the quasi-static sheh is given by and the condition for the magnetosphere insta- 
bihty can thus be rewritten as: 

r ^ 2h-i)(i + ,.^?) „(^,,„,)^,.^ (14) 

TcT 7 cos X 

According to [T7^, when the external gas pressure decreases with radius as P ^ the form of 

the magnetosphere far from the polar cusp can be described to within 10% accuracy as (cos A)"-^^^"^ (here 
A is the polar angle counting from the magnetospheric equator). The instability first appears near the 
equator, where the curvature is minimal. Near the equatorial plane (A = 0), for a poloidal dependence of 
the magnetosphere « (cos A)"'^^ we get for the curvature fcpi?A = 1 + 0.27. The toroidal field curvature 
at the magnetospheric equator is fc^ = 1. The tangent sphere at the equator cannot have a radius 
larger than the inverse poloidal curvature, therefrom kKa = 1-27 at A = 0. This is somewhat larger 
than the value of k,Ra = 7/(2(7 — 1)) = 5/4 = 1.25 ( for 7 = 5/3 in the absence of turbulence or for 
fully isotropic turbulence), but within the accuracy limit0. The contribution from anisotropic turbulence 
decreases the critical temperature; for example, for 7 = 5/3, in the case of strongly anisotropic turbulence 
m|| = 1, m_L = 0, at A = we obtain T jT^r ^ 2, i.e. anisotropic turbulence increases the stability of the 
magnetosphere. So initially the plasma-magnetospheric boundary is stable, and after cooling to T < T^r 
the plasma instability sets in, starting in the equatorial zone, where the curvature of the magnetospheric 
surface is minimal. 

Let us consider the development of the interchange instability when cooling (predominantly Compton 
cooling) is present. The temperature changes as [19], [20] 

dT T-T^ 



(15) 

tc = J-!^3l^ ^ 10.6[s]i?2M-i . (16) 



dt tc 



where the Compton cooling time is 



tc = « IQ.mRlM^,' ■ (17) 

Here rrie is the electron mass, gt is the Thomson cross section, — O.lMc^ is the X-ray luminosity, 
T is the electron temperature (which is equal to the ion temperature since the timescale of electron-ion 
energy exchange here is the shortest possible), is the X-ray temperature and Hm = 0.6 is the molecular 
weight. The photon temperature is = (l/4)Tc„t for a bremsstrahlung spectrum with an exponential 
cut-off at Tcut, typically T^; = 3 — 5 keV. The solution of equation ([T5|) reads: 



T = + {T„ - r,)e-*/*^ . (18) 

We note that Tcr ^ 30keV ^ T^^ ~ 3 keV. It is seen that for t k, 2tc the temperature decreases to T^. 
In the linear approximation the temperature changes as: 

T « T,,(l - t/tc) . (19) 

Plugging this expression into ([T^. we find that the effective gravity acceleration increases linearly with 
time as: 

GM t , , 

fe// ~ ^2-^cosx- (20) 

n-A tc 

Correspondingly, the velocity of matter due to the instability growth increases with time as: 



Ur= I g^ffdt =^*^ COS X. (21) 







^In | 30| . the curvature is calculated to be kR^ Ri 1.34, still within the accuracy limit 
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Here, Unst is the characteristic time of the instabihty which can be expressed in the form: 

The choice of this expression is due to the fact that in the case of rapid coohng, the velocity of matter 
Ur is of the order of the free-fall time m//, and for slow cooling Ur <^ ^ff- We have also defined 
f(u) = Ur/uff < 1, which will be used in the following. Kq is a dimensionless constant of the order of 
unity. 

Plugging tinst into pT]) . we find the velocity obtained by the matter during the time-scale of the 
instability: 

UriUnst)^^ , ^ COSX- (23) 

2 tcj{uY 

Dividing both parts of this equation by u/y and solving for /(u), we get the expression for f{u): 

f{u)-(^y' (^y\cosx)'/'. (24) 



We used here the expression for the free-fall time: 



3/2 



" UffiRA) V2GM' ^ ' 

Then, the characteristic time-scale for the instability can be rewritten in the form: 

i».* = ^^^f^V''(cosx)-V3. (26) 



i^K{RA) \tff 

From this it can be seen that for tc ^ t// , the timescale for the instability is much larger than the 
free-fall time. 



2'^\2Ko)'/' (cosx)-'/' (27) 



iff 

On the other hand, the time-scale of the instability is shorter than the Compton cooling time: 



tc 



2i/2(2ifo)i/3 (^^^ (cosx)-'/' < 1, (28) 



which allows us to use the linear expansion of temperature increase as a function of time time (jl9p . 
The characteristic scale of instability growth is: 

?2 



A = / Urdt = 1 ^ COS X = iurUnst = ^KqRa ■ (29) 

6 Rj^ tc o o 



In this way, during tinst, the scale of the instability becomes comparable to the magnetospheric radius, 
and the settling velocity turns out to be much smaller than free-fall velocity Uff. Clearly, later in the 
non-linear stage of the instability growth the velocity of matter approaches the free-fall velocity. We 
mainly consider the linear stage, since at this stage the temperature is still high enough (although the 
entropy starts decreasing with decreasing radius), and it is in this zone that a toroidal component of the 
magnetic field is formed and effective angular momentum transfer from the magnetosphere to the shell 
can take place. At later stages of instability growth, the loss of entropy is too strong for convection to 
begin. 

Let us estimate the accuracy of our approximation by retaining the second-order terms in the exponent 
expansion. Then the velocity the matter acquires during the instability time tinst is: 



^|f)^^COSx)-/3' 



(30) 
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Clearly, the smaller accretion rate, the smaller the ratio tff/tc, and the better our approximation. 

We note that for the magnetospheric radius in the form ~ cos A" we have tanx = ntanA. Therefore, 
for n ~ 0.27 close to the equator cosx ~ 1 with high accuracy, and we will in the following ignore this 
factor. We also note that in the magnetospheric cusp region cosx — 0, and in this region matter can 
almost not enter the magnetosphere at all. Substituting ([T7)) into (IMl) and then f{u) into definition (|lip . 
we find for 7 = 5/3 the expression for the Alfven radius in this regime: 

Ra « 1.55 X lQ'[cMKl"\l + ^™?)V.(| mO]"'/" (4)'^" ' ^^^^ 

We stress the difference of the obtained expression for the Alfven radius with the standard one, Ra ~ 
^4/7^j^^-2/7^ which is obtained by equating the dynamical pressure of falling gas to the magnetic field 
pressure; this difference comes from the dependence of f{u) on the magnetic moment and mass accretion 
rate in the settling accretion regime. 
The coefficient due to turbulence 

if, = (l + ^m2)V;(|™0 (32) 

is obviously equal to 1 for isotropic turbulence (see the expression for "0 0)j s-nd thus of interest only in 
the case of anisotropic turbulence. 

Plugging (|3T|) into (|24| . we obtain an explicit expression for f{u): 

fin) « 0.39</"iry^^Mf/"M3-o . (33) 

A necessary condition for removal of angular momentum from the magnetosphere via convection is 
the condition of subsonic settling (the Mach number for the settling velocity A4 = Ur/cg < 1), which 
for 7 = 5/3 is reduced to the inequality f{u) < Clearly, this condition is fulfilled for mass 

accretion rates around 10^^ g/s and lower. It is also important to stress that convection in the shell as 
well as removal of angular momentum practically stops working when the mean radial settling velocity 
of the matter Ur becomes higher than the convective velocity u^, i.e. when the convective Mach number 

= Uc/cs ^ rrit becomes smaller than the standard Mach number A4 — Ur/cs- And oppositely, when 
the Mach number of the radial flow becomes smaller than the turbulent Mach number A4 < rrit ~ TOc, 
removal of angular momentum through the shell may take place. When the accretion rate of matter 
through the shell becomes larger than a certain critical value M > the velocity of the accretion flow 
close to the Alfvenic surface may become higher than the sound speed, and a supersonic flow region 
with matter in free fall may form above the magnetosphere. Through this region it is not possible to 
remove any angular momentum from the rotating magnetosphere. In this case, settling accretion is not 
applicable. A Shockwave forms above the magnetosphere and plasma interaction with the magnetosphere 
is described in the scenario studied in e.g. jlOI . Depending on the inhomogeneity of the captured stellar 
wind, the specific angular momentum may be either positive or negative, and thus alternating episodes 
of spin-up and spin-down of the neutron star are possible in the supersonic regime. It is easy to estimate 
the critical X-ray luminosity above which the transition from the subsonic (at low X-ray luminosities) 
to the Bondi-Hoyle-Littleton (at high X-ray luminosities) regime takes place. Indeed, assuming a limit 
for the dimensionless settling velocity of /(u)=0.5 (at which removal of angular momentum through the 
shell is still possible, see further Appendix E), from equation (jSSp . we find the maximum possible value 
of the accretion rate for the settling regime with removal of angular momentum: 

Ml^2K,y'K;'/'^^%\ (34) 

We note that a similar value for the critical accretion rate can be found from a comparison of the Compton 
cooling time to the time-scale for convection close to the Alfven radius. 

To conclude this section, we note that it is not difficult to perform a similar analysis for the velocity 
of matter in the magnetosphere due to radiative cooling of the plasma, for cases when Compton cooling 
is less effective [21 . This scenario may be realized in X-ray pulsars at very low accretion rates when 
the shape of the X-ray beam-pattern changes and the photon beam forms a pencil diagram illuminating 
the magnetospheric cusp. In this way one can explain the episodic jioff-states^^ (with very low X-ray 
luminosity), accompanied with a phase-shift in the X-ray pulse profile [22] as observed in pulsars like e.g. 
Vela X-1. 
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3 Transfer of angular momentum to the magnetosphere 



Let us now consider a quasi-stationary subsonic shell in which accretion proceeds onto the neutron star 
magnetosphere. We stress that in this regime, i.e. the settling regime, the accretion rate onto the 
neutron star is determined by the density at the bottom of the shell (which is directly related to the 
density downstream the bow shock in the gravitational capture region) and the ability of the plasma to 
enter the magnetosphere through the Alfvcnic surface. 

The rotation law in the shell depends on the treatment of the turbulent viscosity (see Appendix B for 
cases when the Prandtl law and isotropic turbulence are applicable) and the possible anisotropy of the 
turbulence due to convection (see Appendix C). In the latter case the anisotropy leads to more powerful 
radial turbulence than perpendicular. In this way, as shown in Appendix B and C, we arrive at a set of 
quasi-power-law solutions for the radial dependence of the angular rotation velocity in a convective shell. 
We shall in the following consider a simple power-law dependence of the angular momentum on radius. 



In Section [S] in applications to real pulsars we will use a quasi-Keplerian law with n = 3/2 as well as 
an iso-angular momentum distribution with n = 2, which in some sense represent limiting cases among 
possible solutions. 

When approaching the bow shock, i? — i?^, and the angular velocity of matter approaches the orbital 
velocity, w — ?> ujb- Close to the bow shock the problem is not spherically symmetric any more since the 
flow becomes very complex (parts of the flow may cause the hot shell to bend, etc.), and the structure 
of the flow can be studied only using numerical simulations. In the absence of such simulations, we shall 
assume that the assumption of an iso-angular momentum distribution is valid up to the front of the bow 
shock located at a distance from the neutron star which we shall take to be the Bondi radius Rb, 



where is the stellar wind velocity at the neutron star orbital distance, and Vorb is the neutron star 
orbital velocity. 

This means that the angular velocity of rotation of matter near the magnetosphere ujm will be related 
to lob via 



(Here the numerical factor J) > 1 takes into account the deviation of the actual rotational law from the 
value obtained by using the assumed power-law dependence near the Alfven radius; see Appendix B and 
C for more detail.) 

Now, let the NS magnetosphere rotate with an angular velocity w* — 2tt/P* where P* is the neutron 
star spin period. The matter at the bottom of the shell rotates with an angular velocity Wm, in general 
difi^erent from uj* . If a;* > ujm, coupling of the plasma with the magnetosphere ensures transfer of angular 
momentum from the magnetosphere to the shell, or from the shell to the magnetosphere if uj* < ujm- 
In the general case, the coupling of matter with the magnetosphere can be moderate or strong. In the 
strong coupling regime the toroidal magnetic field component Bt is proportional to the poloidal field 
component Bp as Bt ^ —Bp{ujm ~ ^*)'t, and \Bt\ can grow to ^ l^pl- This regime can be expected for 
rapidly rotating magnetospheres when uj* is comparable to or even greater than the Keplerian angular 
frequency ljJk{Ra)] in the latter case the propeller regime sets in. In the moderate coupling regime, the 
plasma can enter the magnetosphere due to instabilities on a timescale shorter than the time needed for 
the toroidal field to grow to the value of the poloidal field, so Bt < Bp. 

3.1 The case of strong coupling 

Let us first consider the strong coupling regime. In this regime, powerful large-scale convective motions 
may lead to turbulent magnetic field diffusion accompanied by magnetic field dissipation. This process 
is characterized by the turbulent magnetic field diffusion coefficient rjt . In this case the toroidal magnetic 
field (see e.g. [5] and references therein) is: 



uj{R) - R 



(35) 



RB^2GM/{V^ + vl,,f 




(36) 




uj*)B^ 



(37) 
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The turbulent magnetic diffusion coefficient is related to the kinematic turbulent viscosity as rjt — i^t- 
The latter can be written as: 

i^t =< utk > . (38) 

According to the phenomenological Prandtl law, the average characteristics of a turbulent flow (the 
velocity Ut, the characteristic scale of turbulence U and the shear ujm — w*) are related as: 

■ut - 'tkm - w*| . (39) 

In our case, the turbulent scale must be determined by the largest scale of energy supply to the turbulence 
from the rotation of the non-spherical magnetospheric surface. This scale is determined by the difference 
in velocity between the solidly rotating magnetosphere and the accreting matter that is still not interacting 
with the magnetosphere, i.e. It — Ra, which determines the turn-over velocity of the largest turbulence 
eddies. At smaller scales a turbulent cascade develops. Substituting this scale into equations (|37)) - (l39)) 
above, we find that in the strong coupling regime Bt ~ Bp. 

The momentum of the forces due to plasma-magnetosphere interactions acts on the neutron star and 
changes its spin according to: 

luj* = J ^^^dS = ±K{0)K2-^ (40) 

where / is the neutron star's moment of inertia, w is the distance from the rotational axis and K(9) is 
a numerical coefficient depending on the angle between the rotational and magnetic dipole axes. The 
coefficient K2 appears in the above expression for the same reason as in The positive sign corresponds 
to positive flux of angular momentum to the neutron star (wm > 1^*). The negative sign corresponds to 
negative flux of angular momentum across the magnetosphere (uim < w*). 

At the Alfven radius, the matter couples with the magnetosphere and acquires the angular velocity of 
the neutron star. It then falls onto the neutron star surface and returns the angular momentum acquired 
at Ra back to the neutron star via the magnetic field. As a result of this process, the neutron star spins 
up at a rate determined by the expression: 

loj* = +zMR\uj* (41) 

where z is a numerical coefficient which takes into account the angular momentum of the falling matter. 
If all matter falls from the equatorial equator, z = 1; if matter falls strictly along the spin axis, z — Q. If 
all matter were to fall across the entire magnetospheric surface, then z = 2/3. 

Ultimately, the total torque applied to the neutron star in the strong coupling regime yields 

„2 

Ilo* = ±K{e)K2-^ + zMR\uj* . (42) 
Rj^ 

Using (|lip. we can eliminate M in the above equation to obtain in the spin- up regime (cj^ > lj*) 



Ilo* = 



^A 



V2(7 - 1)(1 + 7m?)^(7, mt)k{e) V Rc 
where — GM / {lo*)'^ is the corotation radius. In the spin-down regime {ujm < lj*) we find 



(43) 



luj* = 



R\ 



^ _ ^ ^lf{u) fRA\ 

^V2(7-l)(l + 7"i?)V'(7,™t)-^W \ 



3/2' 



(44) 



Note that in both cases Ra must be smaller than Rc, otherwise the propeller effect prohibits accretion. In 
the propeller regime Ra > Rc, matter does not fall onto the neutron star, there are no accretion-generated 
X-rays from the neutron star, the shell rapidly cools down and shrinks and the standard lUarionov and 
Sunyaev propeller regime [H], with matter outflow from the magnetosphere, is established. 

During both spin-up and spin-down, the neutron star angular velocity w* almost approaches the 
angular velocity of matter at the magnetospheric boundary, uj* — > ujm{RA)- The difference between 
w* and ujm is small so the second term in the square brackets in (j43l) and (j44p is much smaller than 
unity. Also note that when approaching the propeller regime (Ra — > -Rc), the accretion rate decreases. 
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f{u) 0, the second term in the square brackets vanishes, and the spin evohition is determined solely 
by the spin-down term —K{0)^^/R\. (In the propeller regime, u!m < i^K{RA), < i^*, i^* > i^K{RA) 
). So the neutron star spins down to the Keplerian frequency at the Alfven radius. In this regime, the 
specific angular momentum of the matter that flows in and out from the magnetosphere is, of course, 
conserved. 

Near equilibrium (w* ^ tOm), relatively small fluctuations in M across the shell will lead to very 
strong fluctuations in w* since the toroidal held component can change its sign by changing from +Bp to 
—Bp. If strong coupling actually occurs in nature, this property would be a distinguishing feature of this 
regime. It is known (see eg. [1], [23]) that real X-ray pulsars sometimes exhibit rapid spin- up/spin-down 
transitions not associated with X-ray luminosity changes, which may be evidence that they temporarily 
enter the strong coupling regime. It can not be excluded that the triggering of the strong coupling 
regime may be due to the magnetic field frozen into the accreting plasma that has not yet entered the 
magnetosphere. Accretion of magnetized plasma onto neutron stars is studied in detail in the recent work 

by m- 



3.2 The case of moderate coupling 

The strong coupling regime considered above may be realized in the extreme case where the toroidal 
magnetic field Bt attains a maximum possible value ~ Bp due to magnetic turbulent diffusion. Usu- 
ally, the coupling of matter with the magnetosphere is mediated by different plasma instabilities whose 
characteristic times are too short for substantial toroidal field growth. As discussed above in Section 
12. 1[ the shell is very hot close to the magnetosphere boundary, so without cooling above it the plasma is 
marginally stable with respect to the interchange instability (according to the calculations by ;17 ). 
Let us write down the torque due to magnetic forces applied to the neutron star: 

luj* ^ J ^^dS (45) 

On the other hand, there is a mechanical torque on the magnetosphere from the base of the shell caused 
by the turbulent stresses Wr^: 



WR^mdS, (46) 
where the viscous turbulent stresses can be written as (see the Appendices for more details) 

Wr^ = P^tR^ ■ (47) 
To scpecify the turbulent viscosity coefficient 

Vt = {uck) , (48) 

we assume that the characteristic scale of the turbulence close to the magnetosphere is It ^ Ra, and 
that the characteristic velocity of the turbulent pulsations is determined by the mechanism of turbulence 
in the plasma above the magnetosphere. If there are strong convective motions in the shell, caused by 
heating of its base, then Uc ~ Cg, where Cg is the sound speed. If convection is prohibited, there is 
still turbulence, caused by the shear flow in the shell {u ^ ^/R^, see the Appendices). In this case 
Uc{Ra) ^ u^{Ra) ^ oj*Ra ^ Cs- Obviously, the ratio of the stresses for the different cases turns out 
to be of the order of uj* /ujk{Ra), which for slowly rotating pulsars is around 0.03 — 0.3. Equating the 
torques p5| and (|46)) . we get 

„ duj BtBj, 

We eliminate the density from this expression using the pressure balance at the magnetospheric boundary 
([5]) and the expression for the temperature and make the substitution 

OR ^ -CRT- 

Here we have introduced the dimensionless factor C < 1, characterizing the size of the zone in which 
there is an effective exchange of angular momentum between the magnetosphere and the base of the 
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shell. Then we find the relation between the toroidal and poloidal components of the magnetic field in 
the magnetosphere: 

Bt 7 f Uc \ f oo„^ - uj* \ 

(51) 



Bp ^{^-l)Kt\uff J \C^k{Ra) 



(Here and below we have used the following designations: the free fall velocity Uff = ^ Keple- 

rian frequency at the magnetospheric boundary loj^^Ra) and the correction coefficient due to turbulence 

Kt = (1 + 7m?)V'(7,TOt))- 

Substituting (jSip into (j45p . in case of convection = rricCs (where we have introduced the Mach 
number for convective motions rric), the spin-down rate of the neutron star can be written as: 

/^*=f^)ifX^^^. (52) 
V C / Ra ^k{Ra) 

where Ki is a constant of the order of unity arising from a combination of the parameters in (jSip . In 
this case ([5T|) can be re- written in the form 

Bp-^[C J ujk{Ra)' 

where the geometrical factors arising from the integration of (j45p are included in the coefficient K ^ I. 

If the differential rotation at the base of the shell gives rise to turbulence, Uc ^ = Ld*RA, and the 
expression for spin down takes the form 

-\-)''^wAr-J ^''^ 

where 

«-(^)" 

is the corotational radius (see also [25]). 

Evidently, the breaking torque is in this case smaller by a factor of {Ra/ Rc)^^^ as compared to when 
there are convective motions in the shell. We will call this case the case of weak cuopling. It can easily be 
seen that in this case both the breaking torque and the spin down rate of the neutron star are independent 
of the mass accretion rate (in the limit ujm — we have just Ksd ~ M^/^c; ES)- As will be discussed 
later on, the non-equilibrium pulsar GX 1-1-4 shows during spin-down a negative correlation lu* with 
luminosity |26| . Therefore, we prefer breaking according to ()52p (i.e. with moderate coupling). 

Using the definition of the Alfven radius Ra (fTTj) and the expression for the Kcplcrian frequency lok, 
we can write (I52|l in the form 

IlO* = ZMR\{0Jrn - w*). (56) 

Here the dimensionless coefficient Z is 

^.M^l(pi),,,. ,57) 

f(u) 47 

Substituting in this formula 7 = 5/3 and the expression p4[) . we find 

'7 ~ n Qf-Q /^-^A t'-'i'/ll 2^21/22 ,V-4/ll 1/11 

Z w 0.363 l—jK^' ' M30 • (58) 

Taking into account that the matter that falls onto the neutron star adds the angular momentum 
zMR\ijj* (see Equation (|1T|) above), we get 

luj* = ZMR\{uj,n ~uj*) + zMR\uj* . (59) 

It is obvious that for angular momentum removal from the neutron star through a shell , the coefficient 
Z has to be larger than z. Then the accreting neutron star can episodically spin down (below we will 
explain this statement in more detail). And conversely, ii Z < z, the neutron star can only spin up. 
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If a hot shell is not formed above the magnetosphere (at high X-ray luminosities or low velocity stellar 
winds, see e.g. |27j and references below), then the supersonic or Bondi accretion regime is established 
and no angular momentum can be removed from the neutron star. In this case Z = equation (j59p 
takes the simple form luj* = ZMR\uim, and the neutron star will spin up to a frequency of the order 
of ujx{Ra) regardless of the sign of the difference between the angular momentum of the matter and the 
magnetic field lines uj^ — uj* close to the magnetospheric boundary. Due to conservation of the specific 
angular momentum = ujb{Rb / Ra)'^ ■ Without the presence of a shell the evolution of the angular 
frequency of the neutron star can be described by the equation 

luj* = ZMubR% , (60) 

where the coefficient Z plays the role of the specific angular momentum of the matter. For example, in the 
models of [15] Z ~ 1/4. Numerical modeling of Bondi- Hoyle-Littleton accretion in two-dimensional (e.g. 
[71 HH]) and three-dimensional (e.g. [Hill]) calculations have, however, shown that due to inhomogeneities 
in the stellar wind, accretion becomes non-stationary and the sign of the captured angular momentum 
may change. The sign of Z may thus also be negative and we may observe alternating spin-up and spin- 
down episodes. Such a scenario is often used to explain the observed changes in the sign of the torque 
in accreting X-ray pulsars (see the discussion in [29]). We stress again that this picture is completely 
realistic for X-ray pulsars at high luminosities > 4 x 10"^^ erg/s, when due to the strong Compton cooling 
around the rotating magnetosphere no convective quasi-hydrostatic shell can be formed. 

If a hot shell is indeed formed (at moderate X-ray luminosities less than ~ 4 x 10'^^ erg/s, see ([M]) '). 
the angular momentum from the neutron star can be transferred outside through the convective shell by 
means of turbulent viscosity. Therefore, substituting ujrn from (|36[) and (|59p . we get 

/w* = ZM&lobR^b " " z/Z)MR\uj* . (61) 

This is the main formula that we will use in the following to describe the evolution of the spin of the 
neutron star. 

The dimensionless coefficients in this equation can be calculated using the factor f{u), which is 
included in the expressions for Z and Ra- Thus, the only dimensionless parameter in the model is ^"^^ • 
Below we will show how this coefficient can be determined using observational data from real X-ray 
pulsars. 



4 Spin-up and spin-down of X-ray pulsars 

In this section we will study the dependence of the accelerating and decelerating torques on the accretion 
rate M. We would like to stress, again, that in our case accretion is subsonic and the accretion rate 
is determined by the ability of matter to enter the magnetosphere through the shell. The velocity with 
which the plasma enters the magnetosphere is then mainly dependent on the density at the magnetospheric 
boundary. The density distribution in the shell is on the other hand directly connected to the density of 
matter in the Shockwave region and density variations downstream the shock are thus rapidly translated 
to corresponding variations in the density near the magnetospheric boundary. This means that variations 
of the accretion rate onto neutron stars in binary systems with circular or low-eccentricity orbits should 
be essentially independent of orbital phase, and be mostly determined by variations in the stellar wind. 
In constrast, possible changes in the capture radius Rb (for example due to velocity changes in the 
stellar wind or variations in the orbital velocity of the neutron star) have little effect on the accretion 
rate through the shell, but strongly affect the torques applied to the neutron star (see Equation (pTj) ). 
Equation (|6ip can be rewritten in the form explicitely showing spin-up and spin-down torques: 

lu* = AM^ ~ BM^I^^ . (62) 

For a characteristic value of the accretion rate Mig = Af /lO^^ g/s, the coefficients (not dependent on 
the accretion rate) will be equal to (in CGS units): 

/ ST, \ 2T1 + 3 9 + 6ti 13 

A « 4.22 X 10^i(0.0388)2-"w -r- i^o " -^t M30 
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B « 5.47 X 10^^(1 - z/Z) (^) K,^'''KV''^4>"' (^) " (64) 

(From now on we will assume 7 = 5/3 in all numerical estimates.) The dimensionless factor 5 < 1 takes 
into account the actual location of the gravitational capture radius, which for a cold stellar wind may be 
somewhat smaller than the Bondi radius [31j . The capture radius can also be reduced due to radiative 
heating of the stellar wind by the X-rays from the neutron star (see below). To derive numerical values 
of the coefficients in Equations and we used the expressions for the coefficient Z (|57l) using 

([551) and dni]) for the Alfven radius. 

Below we will study the case Z — z > 0, i.e. S > 0, since in the opposite case only spin-up of the 
neutron star is possible. 



4.1 Equilibrium pulsars 

For equilibrium pulsars we set uj* — and from Equation (j59p we get 

Zaq{ujm- w*) + zuj* =Q . (65) 



Close to equilibrium we may vary (1591) with respect to M . It is convenient to introduce the dimen- 
sionless parameter y = M / Meq , so that close to equilibrium y = 1 . Variations in SM may in general be 
caused by changes in density Sp as well as in velocity of the stellar wind Sv (and thus the Bondi radius) . 
From the continuity equation and taking into account the dependence of f(u) on M in the shell psp . we 
get 

-^^ = ^-3^ (66) 
Let us start by studying variations in the density only. Assuming Rb = const, we find 

'M^^' ^ ^ n"'" + n^^^^-^^A (67) 

Using the expression for Um from ((S5)) and substituting it into (p7)) . we get 

7 '1—^, - ^5^*"^ ^"^l %\y='\ ( Zl^ ^-4/11^6/11 .y-7/11 -12/11 . . 

Now let us keep the density constant and study changes in the velocity only. Then, we have from (|66p 
that bvjv = -(7/33)(5A'//A/. Varying dSg), we get 

5n-3 _ _ 1-1 / '%\y=A f P* \ ^-4/11 6/11/.-7/11 -12/11 . . 

^^^'^ " - ~ ~ [To^ ) [loo's ) ^" ^^is ■ ^^^^ 

A majority of neutron stars in X-ray pulsars rotate close to their equilibrium periods, i.e. on average 
id* = 0. Near equilibrium we get from (1621) in the settling accretion regime: 

- (mM:::!.)^,-'-.-,- {^f . m 

This expression can be reversed to give the equilibrium period for a system if the magnetic field is 
known: 

D ^ 1300[s] X j^2n/ll j^-3n/ll 6n/ll/ Pb \ ,Y-2n/ll/ "sV" ^-71 ^ 

The ratio of pulsar's to Keplerian frequency at the Alfven radius is independent of n and equal to 

'■'* „ r,v 7^3/11 7^-9/22 / P* \ ^ 9/11 ,-,-3/11 /-r,N 

0.27Kg' —— p^'q . (72) 



ujk{Ra) ' " \100t 
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GXX+1 
SXP1062 



y=M/M_, 



Figure 2: An illustration of the dependence of lo* on the dimensionless accretion rate y [(|82p]. In fact 
as y — >■ 0, Cj* approaches some negative value since the neutron star enters the propeller regime at small 
accretion rates. The figure shows the position in the diagram for equilibrium pulsars with y ^ 1 and for 
non-equilibrium pulsars at steady spin-down with y < y^r 



At equilibrium, the ratio between the toroidal and polodial magnetic fields at the Alfven radius (Equation 
([5T|) ) takes the form: 



El\ /i^iA ( z \f Lu* \_ 10fiu)z / UJ* 



(73) 



Substituting f{u) and ([7^ in this expression, we get: 



Bp 



7^10/11 / n 

-0.75.^f4:) .TmT (74) 

eq 



kI"^/^^ VlOOs 



We stress that for slowly rotating accreting pulsars the ratio between the neutron star spin frequency 
and the Keplerian frequency at the Alfven radius is always smaller than unity. Therefore, for typical 
values f{u) ^ 0.3 and z — 2/3 we have Bt/Bp < /u!k{Ra)) < 1, and the pulsars are far from being 

in the propeller regime (see further discussion in Section [6.2p . 

We would like to stress that in the important case n — 2 (iso-angular- momentum distribution), 
the coefficient in the second term in (|68p vanishes, and thus equating Z^q to(|58p we find the value of 
the magnetic moment of the neutron star only from the pulsar equilibrium period and the derivative 
{duj/dy)eq: 

- ^ (Is) i^) i^y '^■r'AT-M-'" ^ (7.) 

For the case n — 2 and a known fi^q we obtain the stellar wind velocity: 

1/4 



^ « 0.53(1 - ./Z.,)-V4^-Vn^3/..^,i/ii^-^3/n (^^J 



(76) 



As will be shown below, for real equilibrium pulsars zjZ^q <C 1, and thus the derived formula gives a 
correct estimate of the stellar wind velocity. Note the weak dependence of the formula on the dimensionless 
constant as well as on the accretion rate. In the framework of our model we may thus, with knowledge of 
the equilibrium spin period P* , the binary period Pf, and with an estimate of the neutron star magnetic 
field determine the stellar wind velocity, without complicated spectroscopic measurements. 

4.2 Non-equilibrium pulsars 

Below we will study the case Z — z > 0, and thus _B > 0, since in the opposite case only spin-up is 
possible. 

First of all, we note that the function cj*(M) reaches a minimum for some Mcr- Differentiating 
equation (I62p with respect to M and equating the achieved expression to zero, we find 

11 

r R 1 
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For M = Mcr the expression lo* reaches an absolut minimum (see. Figl5|). 
It is convenient to introduce the dimensionless parameter 



M , , 

(78) 



Meg 

where M(.q represents the accretion rate at which w* = 0: 

. .ll/2n 

^4, = (79) 

Obviously, 

11 

M.. = Me,(^)^\ (80) 
In other words, uj* has a minimum for a value of the dimensionless parameter of 



l/c.= (^) <1. (81) 



Equation (|62p can be rewritten in the form 



3 + 2ti 



1^* = AMeq'^ y— [l^y-T^j, (82) 
The minimum cj* for y = ycr (i-e. the maximum possible spin-down rate of the pulsar) is 

^^™„ = -Y^Me," y^. (83) 
Now, we vary (|82]) with respect to y: 



= /^(5y) = A^MeFy^'/" (^y^ - l) (^y) . (84) 

Apparently, depending on whether y > ycr or y < ycr, correlated changes of Sld* with X-ray flux should 
have different signs. Indeed, for GX 1+4 in p6] and [32] a positive correlation of the observed SP with 
SAI was found using the CGRO BATSE and Fermi GBM data. This means that there is a negative 
correlation between Suj* and 6M , suggesting y < ycr in this source. 

Let us now consider accreting pulsars in the stage of spin-down (like e.g. GX 1+4 and SXP 1062). 
If the pulsar is spinning down, measurements of the spin-down rate give limits on the parameters of our 
model. From the simple fact that the spin down is stable, using equations ([62]) . ((63)l and ([64]) we may 
obtain a lower limit on the magnetic field in the case of quasi-spherical accretion with uj* < 0, 



/i30>M30,mm~ 0.1(1 -Z/Z) 12 W i^g ' K/ I — - I M^g' 



(85) 



(and thus equation (|70)) is here transformed into an inequality) . We now make use of the fact that during 
spin down there is a maximum possible breaking torque (see equation (j83p ). Inserting the values of the 
coefficients A and B from equations and ([M)) into we find: 

-1.13 X io-"i„d/.i(i - (^) A {^)' {^y {^y . m 



For the accretion rate M — Mcr this expression reaches the numerical value 

U2(1-./Z,W.WA(^)"(|^)' 

(Note the extremely strong dependence on the stellar wind velocity.) 



(87) 
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Then, from the condition \uj*J < luJ^^^^.J 
magnetic field: 



follows a more interesting lower limit on the neutron star 



> fJ-30 



0.94 



sd 



10-i2rad/s 



X -1/2 / X -3/2 



c 



WsJ 



lOOs 



7/8 f -3/8 



Note the weaker dependence of this estimate on the stellar wind velocity as compared to the inequality 

If the accelerating torque can be neglected compared to the breaking torque (corresponding to the 
low X-ray luminosity limit y <^ 1), we find directly from (|52p that for accreting pulsars at spin down, 
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From this we obtain a lower limit on the neutron star magnetic field that does not depend on the 
parameters of the stellar wind nor the binary orbital period: 
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from ((53|l and (|52]) we get: 
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We see from (pij) , that with decreasing M the ratio Bt j Bp increases for reasons well understood — at low 
M the characteristic cooling time for the plasma increases and the toroidal component has time to grow 
to the same strength as the poloidal. Bt can, however, not become larger than Bp due to an instability 
similar to that of a tightly wound spring. Equating Bt = Bp, and using (f9T|) . we find the luminosity 
below which the pulsar enters the strong coupling regime during spin down (see Section 13.11 above) : 
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Below this luminosity in the strong coupling regime the spin-down law becomes i^s^ ^ is^R^^ 
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M6/11: 
(93) 



(Note that when the spin-up torque can be neglected the expression does not contain the - ever so hard 
to determine - velocity of the stellar wind. ) 

For a further decrease of the accretion rate in non-equilibrium pulsars, the Alfven radius will grow 
to the corotation radius and the pulsar may enter a transient state (the propeller regime). From the 
condition co* = ^/GM/R\ we find the accretion rate for this transition: 



M- 



16 
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(94) 



The formulae derived above show that the restrictions on the model become more significant if the 
neutron star magnetic field can be measured independently (for example using spectral cyclotron lines). 
We also would like to stress the fact that measurements of correlated fiuctuations of the spin frequency 
derivative with luminosity during spin down allows us to place the source in a w* — ?/ diagram (see Fig. 
[2]). To the right from the minimum y > j/cr and the correlation positive. To the left y < ycr and the 
correlation is negative. This way we may obtain further limits on the parameters of our model. Below we 
will perform this analysis for the source GX l-f4, in which such correlations where measured [26], |32j . 
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Figure 3: Torque-luminosity correlation in GX 301-2, ili* as a function of BATSE data (20-40 keV pulsed 
flux) near the equilibrium frequency, see |47] . The assumed X-ray flux at equilibrium (in terms of the 
dimensionless parameter y) is also shown by the vertical dotted line. 

5 Application to real X-ray pulsars 

In this Section, as an illustration of the possible applicability of our model to real sources, we will consider 
five particular slowly rotating moderately luminous X-ray pulsars: GX 301-2, Vela X-1, GX l-f4, SXP 
1062 and 4U 2204-1-56. The first two pulsars are close to the equilibrium rotation of the neutron star, 
showing spin- up/spin-down excursions near the equilibrium frequency (apart from the spin- up/spin-down 
jumps, which may be, we think, due to episodic switch-ons of the strong coupling regime when the toroidal 
magnetic field component becomes comparable to the poloidal one, see Section [3TT|) . The third source, GX 
1-1-4, is a typical example of a pulsar displaying long-term spin- up/spin-down episodes. During the last 
30 years, it has shown a steady spin-down with frequency fluctuations (anti-)correlated with luminosity 
(see |32| for a more detailed discussion). Clearly, this pulsar can not be considered to be in equilibrium. 
The pulsar SXP 1062 in the Large Magellanic Cloud as well as the pulsar 4U 2206-1-54 have only been 
observed at steady spin-down. 

5.1 GX 301-2 

GX301-2 (also known as 4U 1223-62) is a high-mass X-ray binary, consisting of a neutron star and an 
early type B optical companion with mass ~ 4OM0 and radius ~ 6Oi?0. The binary period is 41.5 days 
[3 9) . The neutron star is a ~ 680 s X-ray pulsar ,40j , accreting from the strong wind of its companion 
{Mioss ~ lO~^M0/yr, I41j). The photospheric escape velocity of the wind is Vesc ~ 500 km/s. The 
semi-major axis of the binary system is a ~ 17Oi?0 and the orbital eccentricity e w 0.46. The wind 
terminal velocity was found 41 to be about 300 km/s, smaller than the photospheric escape velocity. 

GX 301-2 shows strong short-term pulse period variability, which, as in many other wind- accreting 
pulsars, can be well described by a random walk model [42}. Earlier observations between 1975 and 1984 
showed a period of ~ 700 s while in 1984 the source started to spin up [13]. The almost 10 years of 
spin-up were followed by a reversal of spin in 1993 |44) after which the source has been continuously 
spinning down [45], [46], |47) . Rapid spin- up episodes sometimes appear in the Fermi GBM data on top 
of the long-term spin-down trend [23] . It can not be excluded that these rapid spin-up episodes, as well 
as similar ones observed in BATSE data, reflect a temporary entrance into the strong coupling regime, 
as discussed in Section 2.4.1. Cyclotron line measurements [IS] yield a magnetic field estimate near the 
neutron star surface of Bq w 4.4 x 10"^^ G (fj, = 1/2BqRq = 2.2 x 10^° G cm'^ for the assumed neutron 
star radius i?o = 10 km). 

In Fig. [3] we have plotted uj* as a function of the observed pulsed flux (20-40 keV) according to 
BATSE data (see [47] for more detail). We will consider the neutron star magnetic field in this source to 
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Figure 4: The same as in Fig. [3] for Vela X-1 |55) . 



be known from observations. An estimate of M can be inferred from the X-ray flux provided the distance 
to the source is known, which is generally not the case to a great certainty. We shall assume that near 
equilibrium a hot quasi-spherical shell exists in this pulsar and that the accretion rate is 3 x 10^^ g/s, 
i.e. not higher than the critical value M* ~ 4 x 10^^ g/s [ ([M)) ]. The derivative du)* jdy can be derived 
from the w* - X-ray flux plot, since in the first approximation the accretion rate is proportional to the 
observed pulsed X-ray flux. Near the equilibrium (the torque reversal point with w* =0), we find from 
a linear fit in Fig. ^duj* jdy w 1.5 x 10^^^ rad/s^. 

The obtained parameters (Z, {^^^^ etc.) for this pulsar are listed in Table 1. We note that the toroidal 
component of the magnetic field is much less than the poloidal (the pulsar is far from the strong-coupling 
limit). The stellar wind velocity, determined using the formula (|76l) . is close to the photospheric escape 
velocity. We also note that the value of the parameter describing the coupling between the plasma and 
the magnetosphere i^i/C is of the order of 14, although by its physical sense the coefficient K\ should be 
of the order of 1. This means that the value of the parameter C, which gives the characteristic relative 
size of the region in which transfer of angular momentum from the shell takes place to the magnetosphere 
(or vice versa) has to be of the order of 1/10 (i.e. the characteristic size of the region where angular 
momentum transfer takes place should be approximately 1/10 of the Alfven radius). 

5.2 Vela X-1 

Vela X-1 (=4U 0900-40) is the brightest persistent accretion-powered pulsar in the 20-50 keV energy band 
with an average luminosity of L:r sa 4 x lO^^erg/s [13]. It consists of a massive neutron star (1.88 M©, 
pS] ) and the B0.5Ib super giant HD 77581, which eclipses the neutron star every orbital cycle of ~ 8.964 
days [in]. The neutron star was discovered as an X-ray pulsar with a spin period of ^^283 s [5D], which 
has remained almost constant since the discovery of the source. The optical companion has a mass and 
radius of ~ 23 Mq and ~ 30 Rsun respectively [IH]- The photospheric escape velocity is Vesc ~ 540 km/s. 
The orbital separation is a « 5Oi?0 and the orbital eccentricity e « 0.1. The primary almost fills its 
Roche lobe (as also evidenced by the presence of elliptical variations in the optical light curve, [51]. The 
mass-loss rate from the primary star is 10~^ MqIjx (Nagase et al. 1986) via a fast wind with a terminal 
velocity of ^ 1100 km/s [53], which is typical for this class. Despite the fact that the terminal velocity 
of the wind is rather large, the compactness of the system makes it impossible for the wind to reach this 
velocity before interacting with the neutron star, so the relative velocity of the wind with respect to the 
neutron star is rather low, ~ 700 km/s. 

Cyclotron line measurements [H] yield the magnetic field estimate i?o ~ 3 x 10^^ G (/i = 1.5 x 10^° G 
cm'^ for the assumed neutron star radius 10 km). We shall assume that in this pulsar M ~ 3 x 10^^ g/s 
(again for the existence of a shell to be possible). In Fig. |4]we have plotted Cj* as a function of the 
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observed pulsed flux (20-40 keV) according to BATSE data |55j. As in the case of GX 301-2, from a 
linear fit we find at the spin- up/spin-down transition point dio* /dy ~ 1.2 x 10~^^ rad/s^. 

The obtained parameters for Vela X-1 are listed in Table 1. We note that the velocity of the stellar 
wind as obtained using ((76)) is very close to the observed value of 700 km/s. As in the case of GX 1+4, 
the value of the coupling parameter i^i/C is of the order of 10, i.e. the size of the region for transfer of 
angular momentum between the plasma and the magnetosphere is about 1/10 of the Alfven radius. 

5.3 GX 1+4 

GX 1-1-4 was the first source to be identified as a symbiotic binary containing a neutron star [56 . The 
pulse period is ^ 140 s and the donor is an Mill giant [S^. The system has an orbital period of 1161 days 
[57) . making it the widest known LMXB by at least one order of magnitude. The donor is far from filling 
its Roche lobe and accretion onto the neutron star is by capture of the stellar wind of the companion. 

The system has a very interesting spin history. During the 1970's it was spinning up at the fastest rate 
{Cjsu ~ 3.8 • 10~^^ rad/s) among the known X-ray pulsars at the time (e.g. [43])). After several years of 
non-detections in the early 1980's, it reappeared again, now spinning down at a rate similar in magnitude 
to that of the previous spin-up. At present the source is steadily spinning down with an average spin 
down rate of omega*^ ~ —2.34 x 10^^^ rad/s. The observed spin-reversal has been interpreted in terms 
of a retrograde accretion disc forming in the system [S^ , [59] , [2^ . A detailed spin-down history of the 
source is discussed in the recent paper [32]. Using our model this behavior can, however, be readily 
explained in the framework of quasi-spherical accretion. 

As the pulsar in GX 1-1-4 is not in equilibrium, we use one of the three formulas from Section [4.21 to 
derive a lower limit on the neutron star magnetic field from the observed value of iogd- From ([88)l we get 
t^30 min ~ ^'^{Ki/ C)~^^^ ■ Assuming that the coupling parameter for non-equilibrium pulsars is similar 
to that in equilibrium ones (and thus that the size of the region where transfer of angular momentum 
between the plasma and the magnetosphere takes place is of the order of 1/10 of the Alfven radius, 
C ~ 0.1) we find that n'^o^^i^ ~ 4. 

In this source we also observe anti-correlated variability in spin-down rate versus X-ray luminosity [26j . 
According to the latest Fermi GEM data in the paper [32] it was found that — cj* ^ L'^'^. In our model 
for moderate coupling, Kgd ~ M^^^^, which is very similar to the observed relation. In the earlier BATSE 
observations [26] it was found that —ut* ~ ^2'^^. It can not be excluded that the average luminosity of 
the source was lower at this time. In that case the component Bt could have been closer to Bp, and 
then the expected correlation would have had the form Ksd M^^^^ ~ L^'^'^- Note that in the model 
with weak coupling (with transfer of angular momentum due to turbulence close to the magnetosphere 
[5S]), the breaking torque is less effective by a factor of {Ra/ RcY^"^ and not at all dependent on the 
luminosity. In low-luminosity pulsars the cooling close to the Alfven radius is less effective, which leads 
to the development of convective movements in the shell and the establishment of the moderate coupling 
regime. 

Further, we note that the short-term spin-up episodes, sometimes observed on top of the steady spin- 
down behaviour (at about MJD 49700, see Fig. 2 in [26^ ) are correlated with an enhancement of the 
X-ray ffux, in contrast to the negative frequency-flux correlations discussed above. During these short 
spin-ups, w* is about half the average a;*„ observed during the steady spin-up state of GX 1-1-4 up to 
1980. The X-ray luminosity during these episodic spin-ups is approximately five times larger than the 
mean X-ray luminosity during the steady spin-down. We remind the reader that once M > M, , a free-fall 
gap appears above the magnetosphere, and the neutron star can only spin up. When the X-ray flux drops 
again, the settling accretion regime is re-established and the neutron star resumes its spinning-down. 

5.4 SXP 1062 

This recently discovered young X-ray pulsar in Be/X-ray binary system, located in a supernova remnant in 
the Small Magellanic Cloud. Its rotational period is P* ~ 1062 s and it has a low X-ray luminosity of « 
6 X 10^^ erg/s [60]. The source shows a remarkably high spin-down rate of w* « —1.6 x 10~^^ (rad/s^). Its 
origin is widely discussed in the literature (see e.g. [61], [37]) and a possibly anormously high magnetic 
field of the neutron star has been suggested |63| . In the framework of our model we use more conservative 
limits. Neglecting the spin- up torque (|90|) . we get fx^o > Mao min ~ 20. This shows that the observed spin 
down can be explained by a magnetic field of the order of 10^'^ G, and thus we believe it is premature to 
conclude that the source is an accreting magnetar. 
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Table 1: Parameters for the pulsars discussed in Section [SJ References for the observed pulsar and orbital 
parameters are given in the text as well as values for the wind velocities from measurements of the optical 
components. The parameters Z , Ki/C, and f{u) were derived in Sections 12.31 and [3l Numerical estimates 
are given assuming iso-angular-momentum rotation in the shell {n = 2), moderate coupling between the 
plasma and the magnetic field i5=l,C = l,Cij = l, Kq = 1,7 = 5/3 without turbulence {rat = 0, Kt = 1). 



Pulsars 


equilibrium pulsars 


non- 


equilibrium pulsars 




GX301 - 


- 2 VelaX - 1 


GXl + 4 


SXP1062 


4U2206 + 54 


Measured parameters 


P*{s) 


680 


283 


140 


1062 


5560 




41.5 


8.96 


1161 


- 300t 


19 




300 


700 


200 


- 300* 


350 


M30 


2.7 


1.2 


7 


? 


1.7 


Mi6 


3 


3 


1 


0.6 


0.2 


If |,=i(rad/s^) 


1.5 • 10- 


-12 1.2.10-12 


n/a 


n/a 


n/a 










-2.34- 10-" 


-1.63- 10-" 


-9.4-10-1'' 


Derived parameters 




0.53 


0.57 
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2.6 








Bt/ Bp 


0.17 


0.22 








i?A(cm) 


2-109 1.4-109 








iAj*/uJK{RA) 


0.07 


0.08 








Vw,min (km/ s) 


500 


740 








M30, mm 






/;' « 4 


Mmm ~ 20 


^^'min ~ 3.6 



''' Estimate of the source's position in the Corbet diagram * Estimate of typical wind velocity binary 

pulsars containing Be-stars. 
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5.5 4U 2204+56 



This slowly rotating pulsar has a period of P* = 5560 s and shows a spin-down rate of ujsd ~ —9.4 x 
-[^Q-14 pad/s |64]. The orbital period of the binary system is Pf, ~ 19 days [64], and the measured 
stellar wind velocity is vyy « 350 km/s, abnormally low for an 09.5V j65| optical counterpart. The 
X-ray luminosity of the source is on average L^; ~ 2 x 10'^^ erg/s. A feature in the X-ray spectrum 
sometimes observed around 30 keV can be interpreted as a cyclotron line [ST], [SH], [51]. That gives 
an estimate of the magnetic field of the order oi B ^ (30/11.6) x 1.3 w 3.4 x 10"'^^ G (taking into account 
the gravitational redshift close to the surface 1 + z ^ 1.3), and thus fi^o « 1.7. Using this value of the 
magnetic field and neglecting the accelerating torque, from the formula in (1891) we obtain a lower limit on 
the parameter Ki/ Q > 8, which is very close to the coupling parameter values for the equilibrium pulsars 
Vela X-1 and GX 301-2. If we consider the magnetic field to be unknown (see discussion in 64^), and 
apply the formula like in the case of GX 1+4, assuming moderate coupling with Ki/Q ^ 10, we get 
the limit /i3o > /i3Q ~ 3.6, which is in agreement with standard neutron star magnetic field values. 
Note that using our formulas for equilibrium pulsars would here give a magnetar value for the magnetic 
field [51]. 



6 Discussion 

6.1 Physical conditions inside the shell 

For an accretion shell to be formed around the neutron star magnetosphere it is necessary that the matter 
crossing the bow shock does not cool down too rapidly and thus starts to fall freely. This means that the 
radiation cooling time tcooi must be longer than the characteristic time of plasma motion. 
The plasma is heated up in the strong shock to a temperature 



A,4.i.36xlO^[K](^^^^^^/^ 



Tj,s = — MmT^ « 1.36 X 105[K] ( 77^7^^ ) . (95) 



The radiative cooling time of the plasma is 



3fcT 
2^,„neA 



(96) 



where p is the plasma density, = Yg^p/mp is the electron number density ( fim = 0.6 and Yg « 0.8 for 
fully ionized plasma with solar abundance). A is the cooling function which can be approximated as 

0,T < lO^K 
. l.Ox 10-24^0-55, io4k<T< 105 K 
^^^^>-\ 6.2x 10-i9r-0-6,105K<r<4x lO^K ^^'> 
2.5 X 10-27^0-5, T > 4 X lO^K 

m, m- 

Compton cooling becomes effective from the radius where the gas temperature T, determined by the 
hydrostatic formula is lower than the X-ray Compton temperature T^. The Compton cooling time 
(see (dZJ)) is: 

« 1060[s]M,V . (98) 

Above the radius where Tx = T, Compton heating dominates. Taking the actual temperature close to 
the adiabatic one [(jl])], we find ^ 2 x 10^° cm. We note that both the Compton and photoionization 
heating processes are controlled by the photoionization parameter ^ [72] , [73] 

(99) 

In most part of the accretion flux, n ~ R~^/^, so ^ ^ R~^/^ and independent of the X-ray luminosity 
through the mass continuity equation. We derive a characteristic value for ^: 

C « 5 X lQ''f{u)R-^'^ . (100) 
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If Compton processes were effective everywliere, this high value of the parameter ^ would imply that the 
plasma is Compton-heated up to keV-temperatures out to very large distances ~ 10^^ cm. However, at 
large distances the Compton heating time becomes longer than the characteristic time of gas accretion: 

= ^^^^^ - 20/(.)A/ra^i?^ (101) 

^accr 

which shows that Compton heating is ineffective. The gas temperature is determined by photoionization 
heating only and the gas can only be heated up to T^ax ~ 5 x 10^ K [72], which is substantially lower 
than 1^-3 keV. 

The effective gravitational capture radius corresponding to the sound velocity of the gas in the 
photoionization-heated zone is 

Everywhere up to the bow shock photoionization keeps the temperature at a value ~ Tmax- The sound 
velocity corresponding to Tmax is approximately 80 km/s. If the stellar wind velocity exceeds 80 km/s, 
a standard bow shock is formed at the Bondi radius with a post-shock temperature given by (1951) . If the 
stellar wind velocity is lower than this value, the shock disappears and quasi-spherical accretion occurs 
from R*g. The photoionization heating time at the effective Bondi radius 3 x 10^^ cm is 



(3/2)fcr,„ax/A* 

{hVeff - Ceff)njaeffC 



« ^^^ l'^/^;'^^^max/A^r» ^ ^ ^ 10%]M^I . (103) 



(here hveff ~ 10 keV is the characteristic photon energy, C is the effective photoionization potential, 
(Te// ~ lO"^** cm^ is the typical photoionization cross-section and = L/(47ri?^/iz/e//c) is the photon 
number density). The photoionization to accretion time ratio at the effective Bondi radius is then 



« 0.07/(zi)MiV . (104) 

'^accr 

At wind velocities Vw > 80 km/s the bow shock stands at the classical Bondi radius Rb inside the 
effective Bondi radius i?^ determined by (|102[) . The cooling time of the shocked plasma at i?^ expressed 
through the wind velocity is: 

tcooi - 4.7 X 10%]M-^\^^-^ . (105) 

The photoionization heating time in the post-shock region can also be expressed through the stellar wind 
velocity: 

tp^ w 3.5 X 10^[?]Aff6^wf . (106) 

A comparison of these two characteristic timescales implies that for low wind velocities radiative cooling 
becomes important and the source enters the regime of free-fall accretion with conservation of specific 
angular momentum. 

Thus, for low wind velocities the plasma behind the shock cools down and starts to fall freely. As the 
cold plasma approaches the gravitating center, photoionization heating becomes important and rapidly 
heats up the plasma to Tmax ~ 5 x 10^ K. Should this occur at a radius where Tmax < GM/{TIR), 
the plasma continues its free fall down to the magnetosphere, still with the temperature Tmax, with the 
subsequent formation of a shock above the magnetosphere. However, if Tmax > GM/ {TZR), settling 
accretion will work even for low wind velocities. 

For high-wind stellar velocities > 100 km/s, the post-shock temperature is higher than Tmax, 
photoionization is unimportant, and the settling accretion regime is established if the radiation cooling 
time is longer than the accretion time. From a comparison of these timescales, we find the critical 
accretion rate as a function of of the wind velocity below which the settling accretion regime is possible: 

M|g < O.Uvf^ . (107) 

Here we stress the difference of the critical acccretion rate AI^ and Mt, derived earlier. For M > A'fi 
he plasma rapidly cools down in the gravitational capture region and free-fall accretion begins (unless 
photoionization heats up the plasma above the adiabatic value at some radius), while at M > ~ 
4x 10^^ g/s, determined by (p4|) a free-fall gap appears immediately above the neutron star magnetosphere. 
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6.2 On the possibility of the propeller regime 

The very slow rotation of the neutron stars in X-ray pulsars considered here (GX 1+4, GX 301-2, Vela 
X-1, SXP 1062, 4U 2204+56) with u:*{Ra) < i^K{RA) (see Table 1) makes it hard for these sources to 
enter the propeller regime where matter is ejected with parabolic velocities from the magnetosphere and 
the neutron star spins down. 

Let us therefore start with estimating the important ratio of viscous tensions (~ BtBp) to the gas 
pressure (~ B^) at the magnetospheric boundary. This ratio is proportional to Bt/Bp (see ([75)1 ) and is 
always much smaller than 1 (see Table 1), i.e. only large-scale convective motions where the characteristic 
hierarchy of eddies scales with radius can be established in the shell. 

When uj* > lok{Ra), a centrifugal barrier is formed and accretion stops (the propeller regime). In 
that case the maximum possible braking torque is —^^/R\ due to the strong coupling between the 
plasma and the magnetic field. Note that in the propeller state, interaction of the plasma with the 
magnetic field is by strong coupling, i.e. the toroidal magnetic field component Bt is comparable to 
the poloidal one Bp. It can not be excluded that a hot iso-angular- momentum envelope could exist in 
this case as well, which would then remove angular momentum from the rotating magnetosphere. If the 
characteristic cooling time of the gas in the envelope is short in comparison to the falling time of matter, 
the shell disappears and one can expect the formation of a 'storaging' thin Keplerian disc around the 
neutron star magnetosphere [!74 . There is no accretion of matter through such a disc. It only serves to 
remove angular momentum from the magnetosphere. 

6.3 Effects of the hot shell on the X-ray energy and power spectrum 

The spectra of X-ray pulsars are dominated by emission generated in the accretion column. The hot 
optically thin shell produces its own thermal emission, but even if all gravitational energy were released 
in the shell, the ratio of the X-ray luminosity from the shell to that of the accretion column would be 
about the ratio of the magnetospheric radius to that of the neutron star, i.e. one percent or less. In 
reality, the luminosity from the shell is much smaller. The shell should scatter X-ray radiation from the 
accretion column, but for this effect to be substantial, the Comptonization parameter y must be of the 
order of one. The Thomson depth in the shell is, however, very small. Indeed, from the mass continuity 
equation and ((3T|) for the Alfven radius and (|33|) for the factor f{u), we get: 

i-Rb 

TT= ne{R)(TTdR w 3.2 x lO^^A/fg^^Mao • 
JRa 

Therefore, for the characteristic temperatures near the magnetosphere (see [(|4])]) the parameter y is 

ikT „ 
y = -TT ~ 2.4 X 10-=^ . 

This means that the X-ray spectrum, formed in the region of energy conversion close to the surface of 
the neutron star is not expected to be significantly altered by scattering in the hot shell. 

Large-scale convective motions in the shell introduce an intrinsic time-scale of the order of the free- 
fall time that could give rise to features (e.g. QPOs) in the power spectrum of variability. QPOs were 
reported in some X-ray pulsars (see [75] and references therein). However, the expected frequencies of 
any QPOs arising in our model would be of the order of mHz, much higher than those reported. 

A stronger effect could be the appearance of a dynamical instability in the shell due to increased 
Compton cooling and hence increased mass accretion rate through the shell. This may lead to a complete 
collapse of the shell triggering an X-ray outburst with duration similar to the free-fall time scale of the 
shell (~ 1000 s). Such transient behaviour is observed in the supergiant fast X-ray transients (SFXTs) 
see [7S].) The possible development of such a scenario depends on the specific parameters of the shell 
and needs to be further investigated. 

6.4 Can accretion discs (prograde or retrograde) be present in these pulsars? 

Our analysis of the sample of pulsars in Section [5] suggested that in a convective shell an iso-angular- 
momentum distribution is the most plausible. Therefore, we shall below consider only this case, i.e. using 
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the rotation law uj ^ R ^. As follows from (|6ip. at w* = the equilibrium angular frequency of the 
neutron star is 

We stress that such an equilibrium in our model is possible only when a shell is present. At high accretion 
rates M > M» ~ 4 x 10^^ g/s accretion proceeds in the free-fall regime (with no shell present). 

The equilibrium period for an X-ray pulsar in the quasi-spherical settling accretion regime can be 
derived using the formula dTT]): Peq ^ 1300[s]^3o/"(Pf,/10d)Mi"g*/"D| . 

For standard disc accretion, the equilibrium period is 

Pe^M-l^fili'M-,'^', (109) 

and the long periods observed in some X-ray pulsars can thus, if a is disc present, be explained only 
assuming a very high magnetic field of the neutron star. Retrograde accretion discs are also discussed 
in the literature (see e.g. [29] and references therein). Torque reversals produced by temporary forming 
retrograde discs can in principle lead to very long periods for X-ray pulsars even with standard magnetic 
fields. Such retrograde discs could be formed as a result of inhomogeneities in the captured stellar wind 
[8l[9]. The scenario could, in principle, work for pulsars at high accretion rate, too high for a hot envelop 
to form. 

In the case of GX 1-1-4, however, it is highly unlikely to observe a retrograde disk on a time scale 
much longer than the orbital period (see a more detailed discussion of this issue in [32]). For both GX 
301-2 and Vela X-1, the observed positive torque-luminosity correlation (see Figs. |3]and|31) also rules out 
a retrograde disc in any of these systems. 

To conclude this discussion section, we should mention that in reality, all pulsars (including those 
considered here) demonstrate a complex quasi-stationary behaviour with dips, outbursts, etc. These con- 
siderations are beyond the scope of this paper and definitely deserve further observational and theoretical 
studies. 



7 Conclusions 

In [T^ we presented a theoretical model for quasi-spherical subsonic accretion onto slowly rotating mag- 
netized neutron stars. In this model the accreting matter is gravitationally captured from the stellar 
wind of the optical companion and subsonically settles down onto the rotating magnetosphere forming 
an extended quasi-static shell. This shell mediates the angular momentum removal from the rotating 
neutron star magnetosphere by large-scale convective motions. Depending on the angular velocity of the 
rotating matter close to the magnetospheric boundary this type of accretion can cause the neutron star 
to either spin up or spin down. 

A detailed analysis and comparison with observations of the two X-ray pulsars GX 301-2 and Vela X-1, 
both demonstrating positive torque-luminosity correlations near the equilibrium neutron star spin period, 
shows that the convective motions are most likely strongly anisotropic, and the rotational velocities in 
the shell lu ~ have a near iso-angular-momentum distribution. We note that a statistical analysis of 
long-period X-ray pulsars with Be-components in SMC by [33] also favored the rotation law w ~ R^^. 

The accretion rate through the shell is determined by the ability of the plasma to enter the magneto- 
sphere. The settling regime of accretion which allows angular momentum removal from the neutron star 
magnetosphere can be realized for moderate accretion rates M < A/* ~ 4 x 10^^ g/s. At higher accretion 
rates a free-fall gap above the neutron star magnetosphere appears due to rapid Compton cooling, and 
accretion becomes highly non-stationary. 

From observations of the spin- up/spin-down rates (the angular rotation frequency derivative uj* , or 
duj* /DM near the torque reversal) of slowly rotating equilibrium X-ray pulsars with known orbital periods 
it is possible to determine the main dimensionless parameters of the model, as well as to estimate the 
magnetic field of the neutron star. Such an analysis revealed a good agreement between magnetic field 
estimates obtained using our model and those derived from cyclotron line measurements for the pulsars 
GX 301-2 and Vela X-1. 

Using measurements of the spin period and the orbital period together with an estimate of the neutron 
star magnetic field fx, our model furthermore offers a possibility to estimate the stellar wind velocity of 
the companion, without the need for complicated spectroscopic measurements. 
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For non-equilibrium pulsars there is a maximum possible spin-down rate, depending on the spin period 
P* , the orbital period Pb, the neutron star magnetic field /i and the wind velocity v^. For such pulsars it 
is possible to estimate a lower limit on the neutron star magnetic field using the observed spin-down rate 
and X-ray luminosity. For the pulsars GX 1-1-4, SXP 1062, 4U 2206-1-54 investigated here, our estimates 
are all in agreement with standard field values and observed cyclotron line measurements. 

In our model for quasi-spherical subsonic accretion, long-term spin- up/spin-down as observed in some 
X-ray pulsars can be quantitatively explained by a change in the mean mass accretion rate onto the 
neutron star (and the corresponding mean X-ray luminosity). Clearly, these changes are related to the 
stellar wind properties. 

The model also predicts the specific behaviour of the variations in Suj*, observed on top of a steady 
spin- up or spin-down, as a function of mass accretion rate fiuctuations SM . There is a critical accretion 
rate Mcr below which an anti-correlation of Su* with 6M should occur (the case of GX 1-1-4 at the steady 
spin-down state currently observed), and above which Suj* should correlate with 6M fiuctuations (the 
case of Vela X-1, GX 301-2, and GX 1-1-4 in the steady spin-up state). The model explains quantitatively 
the relative amplitude and the sign of the observed frequency fiuctuations in GX 1-1-4. 



A The structure of a quasi-spherical rotating shell with accre- 
tion 

A.l Basic equations 

Let us first write down the Navier-Stokes equations in spherical coordinates i?, 9, (f>. Due to the huge 
Reynolds numbers in the shell (~ 10^^ — 10^® for a typical accretion rate of 10^^ g/s and magnetospheric 
radius ^ 10^ cm), there must be strong turbulence. In this case the Navier-Stokes equations are usually 
called the Reynolds equations. In the general case, the turbulent viscosity may depend on the coordinates, 
so the equations take the form: 

1. Mass continuity equation: 

dp 1 9 . o X 1 d , . ^ . 1 dpuA, ^ / A 1 X 

T^+min^ P^r) + iT^T^Tj {sine pue) + „ . „ ^ ,^ - 0. A.l 
ot R'^ oR ^ ' R sm oO R sm 0(p 

2. The i?-component of the momentum equation: 

dUr ^ dUr , Ue dUr , U0 dur + GM 

3. The 0-component of the momentum equation: 

dug dug Ugdug dug UrUg - U^^COt 9 

^ + + + R ^""^ (^-'^ 

4. The 0-component of the momentum equation: 

du^ du^ ugdu^ du^ u^u^ + u^ugcoi9 

-dt + ai? + i? ^ + R^9^ + R = ^^-^^ 

Here the force components (including viscous force and gas pressure gradients) read: 

^-^ -hl (^^^-) - ^rI9 - ^1^- - ^ - ^ (^•^) 

- (^'^-) + ^^9 + ^1^^^ - (^-^^ 

- i^l? (^^^^-) + ^^9 + ^^^^^ (^-^^ 

The components of the stress tensor include a contribution from both the gas pressure Pg (assumed to 
be isotropic) and the turbulent pressure P* (generally anisotropic). In their definition we shall follow the 
classical treatment by Landau and Lifshitz |77j but with the inclusion of anisotropic turbulent pressure: 
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f)n, 2 

Wrr = -P, - P|j^ + 2pvt - -pi^tdWu (A.8) 

Wee = -P, Pie + ^P-t + ^) - ^P-^divu (A.9) 

W,, = -P, - PU + 2p.. + p + ^ j - 3^^*^^^" ^^-lO) 

In our problem the anistropy of the turbulence is such that Prr ~ P|| , P^e = P*^^ — P± ■ The turbulent 
pressure components can be expressed through turbulent Mach numbers and will be given in Appendix 
E. 

divu in spherical coordinates is; 

divu = -L^ (R\^) + -1^4 + ir^^- (^-14) 

R^oR^ ' RsmO oO RsmO ocp 

A. 2 Symmetries of the problem 

d . d 

We shall consider axially-symmetric = 0), stationary (— = 0), and only radial accretion (ug = 0). 
Under these conditions, from the continuity equation (jA.ip we obtain: 

M = A-nR^puR = const . (A. 15) 

The constant here is determined from the condition of plasma leakage through the magnetosphere. 

Let us rewrite the Reynolds equations under the above assumptions. The i?-component of the mo- 
mentum (jA.2p equation becomes: 

/ duR ul\ GM 1 d . 2 ^ , 1 9 ,,,, . Wee W^^ 

P [^^^ - P J - 'Pl^ + R-^dR + [WRss^ne) - ^ - ^ (A.16) 

The ^-component of the momentum equation; 

u^cot^ 1 d , n s 1 d Wen 

The 0-component of the momentum equation; 

The components of the stress tensor with anisotropic turbulence take the form; 

W„,^-P,-Pi-P..{-^-'^) (A.19) 

W..^-P,-Pi + l,..(±-!l^) ,A,20) 
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(A.22) 



The main problem in describing gas dynamic flows with turbulence is in finding the kinematic viscosity 
parameter vt ■ As is well known, in the case of laminar flows the viscosity parameter v is dependent only 
on the properties of the medium (liquid or gas). When turbulence is present, however, this coefflcient 
is determined also by the macroscopic properties of the flow itself. The are some empirical relations 
which in principle can be verified experimentally. Most often the so called turbulent mixing length It 
is introduced. Furthermore L. Prandtl in his works introduced for plane-parallel shear flows (along the 
X-axis to be specific) the relation between the turbulent mixing length It, the velocity of the turbulent 
flow Ut and the characteristic amount of shear in the direction perpendicular to the average flow {z): 



du 
dz 



(A.25) 



where Co ~ 1 is a universal dimensionless constant, the exact numerical value of which should be deter- 
mined from a theory that currently does not exist. In this way, the dependence of the turbulent stresses 
on the shear value becomes quadratic: 

du\^ 



and a non-linearity is formed which in the general case makes the problem a lot more difficult. 

First, we consider the general Prandtl law for turbulent viscosity in the case of an axi-symmetric flow. 
In the case of strong anisotropy there is one more empirical law for describing the turbulent viscosity, the 
Wasiutyhski-law (see below), which is not reduced to the Prandtl law in the case of isotropic turbulence. 
This more general case for anisotropic turbulence will be discussed separately in Appendix C. 



B Structure of the shell in the case of turbulent viscosity ac- 
cording to the Prandtl law 

B.l The empirical Prandtl law for ctxisymmetric flows with turbulent viscos- 
ity 

. We consider an axisymmetric flow with a very large Reynolds number. By generalizing the Prandtl 
law for the turbulent velocity obtained for plane parallel flows, the turbulent velocity scales as Ut ~ 
ltR{dui/dR). From the similarity laws of gas-dynamics we assume k ~ R, so 



Ut = CiR^ 



dcu 



OR 



(B.l) 



We note that in our case the turbulent velocity is determined by convection, and thus Ut ^ 0.5u// 
(see Appendix D). This implies that the constant Ci scales as 



Ci ~ ut/{u^), 



(B.2) 



and can be very large since (it^) <C ut- 

The turbulent viscosity coefficient thus reads: 



Ut = {uth) = 



dtu 
dR 



(B.3) 
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Here C2 ~ 1/3 is a numerical factor originating from statistical averaging. Below we shall combine Ci 
and C2 into the new coefficient C — C1C2, which can be much larger than unity. 
For such a viscosity prescription the turbulent stresses M^jj^ are equal to 

B.2 The angular momentum transport equation 

A similar problem (that of a rotating sphere in a viscous fluid) was solved in Landau and Lifshitz |77| . 
They showed that the variables here become separated and u^{R, 9) = u^{R) sin 6*. Note that the angular 
velocity lo{R) = u^{R)/R is independent of the polar angle 6. Our problem is different from that of the 
sphere in a viscous fluid in several respects: 1) there is a force of gravity present, 2) the turbulent viscosity 
varies with R and can in principle depend on 9, and 3) there is radial motion of matter (accretion). These 
differences lead, as will be shown below, to the radial dependence u^{R) cx _R~^/^. (We recall that for a 
rotating sphere in a viscous fluid oc R^^). 

Let us start with solving (IA.18|) . First, we note that for u^{9) ^ sinfl, according to (jA.23|) . Wgrj, = 0. 
Further, making use of the continuity equation (|E.8|) and the definition of angular velocity, we rewrite 
()A.18|) in the form of angular momentum transfer by viscous forces: 

We rewrite equation (jA.24|) using the derivative of the angular velocity: 

WR4 = pvtR^sm9. (B.6) 

Substituting this expression into (|B.5|) and integrating over i?, we get 

MljR^ ^AirpiytR'^^ + D, (B.7) 
oR 

where D is an integration constant. This equation for angular mometum transport by turbulent viscosity 
is similar to that of disc accretion 2 , but different due to spherical symmetry of our problem. 

The left part of (|B.7I) is simply advection of specific angular momentum averaged over the sphere 
(1/2 ujR^ sin^ 9sm9d9 = I/SluR"^) by the average motion toward the gravitational center (accretion). 
M is negative as well as The first term on the right describes transport of angular momentum 

outwards by turbulent viscous forces. 

The constant D is determined from the equation 

(B.8) 

V C / Ra ^k{Ra) 

(see (15^ in the text). We consider accretion onto a magnetized neutron star. When D < 0, the advection 
term in the left part of (jB.7l) dominates over viscous angular momentum transfer outwards. Oppositely, 
when Z? > 0, the viscous term in the right part of (IB.7|) dominates. In the case oi AI — Q (no plasma 
enters the magnetosphere) , there is only angular momentum transport outwards by viscous forces. 
Now let us rewrite (IB.8|) in the form 

D-{'^]k,^R\^-^ (B.9) 



and use the pressure balance condition 



PiRA) = P,{Ra){1 + Wt) = -^C- (B.IO) 



Stt 2w R\ 



Using the mass conitnuity equation in the form 



\M\ = AnR'^pfiu)^/GM/R, 
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and the expression for the gas pressure ([5]), we write the integration constant D/\M\ in the form 

— ^ = — ^ w 7,TOt) ^ (l + 7mf). (o.ll) 

Let us consider the case when the neutron star rotates close to equihbrium Co* = 0. In this case 
according to ([55)) 

Wm-a;* = -— w*, (B.12) 
and thus using definition of Z [ (j57p ]. we obtain: 

D 

W\ 



zR\uj* . (B.13) 



We would like to stress that the value of the constant D is fully determined by the dimensionless specific 
angular momentum of matter at the Alfven radius z. 

B.3 The angular rotation law inside the shell 

Let us now use (jB.7p to find the rotation law uj{R). At large distances R ^ Ra (we would like to remind 
the reader that Ra is the bottom radius of the shell), the constant D is small relative to the other terms, 
so we can set D « 0. Thus, to obtain the rotation law we shall neglect this constant in the right part of 
(|B.7p . Next, we substitute (|B.3I) and make use of the solution for the density (which, as we shall show 
below, remains the same as in the hydrostatic solution) 

p{R)^p{R^)(^^^"^ (B.14) 

in equation (|B.7p to obtain: 



M 

After integrating this equation, we find 



uR^^A^p{Ra)\^-^\ CR'[^] . (B.15) 



where 

\ML_^ (B.ir) 

Atip{Ra)CRY 

and Di is some integration constant. In (IB.16P we use only the positive solution (the minus sign with 
constant Di > would correspond to a solution with the angular velocity growing outwards, which is 
possible if the pulsar has a very long spin period, i.e. almost does not rotate at all). If Di ^ 0, at large 
R 3> Ra (in the zone close to the bow shock) the solid body rotation law would lead to w — const w^. 
(However, we remind the reader that our discussion is not applicable close to the bow shock region.) At 
small distances from the Alfvenic surface the effect of this constant is small and we shall neglect it in the 
calculations below. Then we find 

.jR)^'- l^^l J^T (B.18) 
^ ' 9 ^'kp{Ra)CR\ \ R J ^ ' 

i.e. the quasi-Keplerian law a;(i?) = uJm{RA/R)^^^ ■ The value ujm in the solution given by (jB.18[) is 
obtained after substituting AI from the continuity equation at i? = Ra into (jB.18[) : 

c„,^Cju;{Ra)^^Cj^-^^^^. (B.19) 
y Li Ha 

(Here we have introduced the correction factor w > 1 to account for the deviation of the exact solution 
from the Keplerian law close to Ra)- 
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As uji(Ra) is smaller than the free-fall velocity, the above formula implies that LOm < ujk{Ra), 
lower than the Keplerian angular frequency. For self-consistency the coefficient C in the Prandtl law is 
determined, according to (jB.19l) , by the ratio of the radial velocity to the rotational velocity of matter 

^^4^MM.!^M^. (B.20) 
9 LJ,nRA 9 u^{Ra) 

We note that this ratio is independent of the radius R and is actually constant across the shell. Indeed, 
the radial dependence of the velocity un follows from the continuity equation with account for the density 
distribution (|B.14|) 

Ur(,R)=Ur{RA)\-^\ . (B.21) 

For a quasi-Keplerian law u^{R) ^ so the ratio Ur/u^ is constant. 

Finally, the angular frequency of the rotation of the shell near the magnetosphere Wm is related to 
the angular frequency of the motion of matter near the bow-shock as 

Rs^'" 



w„ = ULOB y-^j . (B.22) 

In fact, when approaching Ra, the integration constant D (which we neglected at large distances R 3> Ra) 
should be taken into account. The rotational law will thus somewhat differ from a quasi-Keplerian close 
to the magnetosphere. 

We stress the principal difference between this regime of accretion and disc accretion. For disc 
accretion the radial velocity is much smaller than the turbulent velocity, and the tangential velocity 
is almost Keplerian and is much larger than the turbulent velocity. The radial velocity in the quasi- 
spherical case is not determined by the rate of the angular momentum removal. It is determined only 
by the "permeability" of the neutron star magnetosphere for infalling matter. In our case we assume 
that the radial velocity is of the order of the velocity of convective motions in the shell. The tangential 
velocity for the obtained quasi-Keplerian law is much smaller than the velocity of the convective motions. 
Note also that in the case of disc accretion the turbulence can be parametrized by only one dimensionless 
parameter a ~ ul/uj. with < a < 1 [2]. The matter in an accretion disc rotates differentially with 
a supersonic (almost Keplerian) velocity, while in our case the shell rotates differentially with a clearly 
subsonic velocity at any radius, and the turbulence in the shell is essentially subsonic. Also, our case with 
an extended shell is of course strongly different from the regime of freely falling matter with a standing 
shock above the magnetosphere [T7] . 

B.4 The case without accretion 

Now let us consider the case where the plasma can not enter the magnetosphere and no accretion onto 
the neutron star occurs. This case is similar to the subsonic propeller regime considered by Davies and 
Pringle [TT]. (|B.7[) then takes the form: 

Q = ATipvtR^^ + D. (B.23) 
oR 

(Remember that the constant D is determined by the spin-down rate of the neutron star, D — < 0). 
Solving this equation as above, we find for the rotation law without accretion: 



Ra\ 
R ) 

where 



7/4 

uj{R) = w„ ( ^ ) , (B.24) 



77tp{RaMRa)R\ 

From (lR3l) we find 



(B.25) 



MRA) = lCiOmR\: (B.26) 
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and thus for Wm we obtain: 



I\uj* 



1/2 



However, Lo,n is also related to the bow-shock region parameters as 



Ra 



(B.27) 



(B.28) 



which can in principle be used to further study this case, which we shall not do here. 



C Structure of the shell and angular rotation law in case of 
turbulent viscosity according to Wasiutynski 

Prandtl's law for viscosity that was used above relates the scale and velocity of turbulent pulsations with 
the average angular velocity and is commonly used when the turbulence is generated by the shear itself. 
In our problem, the turbulence is initiated by large-scale convective motions in the gravitational field. 
Due to convection, strong anisotropic turbulent motions may appear (the radial dispersion of chaotic 
motions could be much larger than the dispersion in the tangential direction), and Prandtl's law may 
thus be inapplicable. 

Anisotropic turbulence is much more complicated and remains poorly studied. As a first step, we 
may adopt the empirical law for Wucf) as suggested by Wasiutynski [75] : 



Wr^ = 2p{-Vt + Vr)^^ + VrpR 



dco 
dR' 



(C.l) 



where the radial and tangential kinematic viscosity coefficients are 

iyr^C\i{\ul\)R 
j^t^C^{\u'J)R 

respectively. The dimensionless constants C|| and C± are of the order of one. In the isotropic case = i^t, 
^ doj/dR, and in the strongly anisotropic case Vr ^ i^t, Wr^ ~ d{ujR^)/dR. Using these definitions, 
let us substitute (IC.ip into (jB.sp . and after integration over R rewrite the latter in the form: 



cjR^ I 1 



2C±(l^^il) 



\Ur\ 



Rd{ujR^ 



dR 



D 

W\ 



{C.2) 



We note that due to self-similarity in the shell structure 



R and thus the ratios 



and {\u^^\)/uj, are constant. In this case the obvious solution to the above equation reads: 



ujR^ 



D 



1 



\M\ 1 - 2C 



LOBRi 



D 



1 



|M| 1 _ 2C 



-Kb \ C|i<i"i|i> 
R 



1~2C, 



(C.3) 



(here the integration constant is defined as such that uj{Rb) = lob)- 

Now let us consider the equilibrium situation where lo* = 0. In this case, as we remember. 



D 

W\ 



~zuj*R\ ,uJm 



(1 - z/Z)u* 



First, let us consider the case of strongly anisotropic, almost radial turbulence where {\i 



this case, the specific angular momentum at the Alfven radius is 



LOmR, 



Rb\<=^\\<H^> 



I- z/Z \ \Ra 



C||(|U,|> 



= 0. In 



(C.4) 
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From this we see that in the case of very weak accretion (or, in the hmit, when there is no accretion 
through the magnetosphere at ah), <^ C||(|w|||), and an almost iso-angular-momentum distribution 
of rotational velocities in the shell is formed. 

The next case is where the amount of anisotropy is such that Cj^(|u*Ll)/|ur| = 1/2. Then we have a 
strict iso-angular-momentum distribution in the shell: — uibR%- 

If the turbulence is fully isotropic C±{\u*^\) = C||(|u||) = C{\u*\). Denoting e = |ur|/(C'(|u*|)), we 
find: 



1 



1 \ (Ra 



1-z/Z J V2/e- 1 / \ \Rb 



ubRI ( ) . (C.5) 



Note that if e — > (and there is no accretion through the magnetosphere), LOm — > ws, and we have 
solid-body rotation without accretion (cf. the first case above!). For e = 3/2, a near quasi-Keplerian 
angular rotation distribution may be established. We remind the reader that a similar quasi-Keplerian 
distribution was obtained in Appendix B above with the use of the Prandtl law for isotropic turbulent 
viscosity. In that case, this was the only solution. Here, in contrast, the quasi-Keplerian law is only one 
particular case of the general solution obtained using Wasiutyhsky's prescription for anisotropic turbulent 
viscosity. 

As we have shown in the main text, a quasi-Keplerian rotation law is not favored by observations. 
We therefore conclude that the most likely velocity distribution in the shell is the near iso-angular- 
momentum one with anisotropic turbulence initiated by convection. Note that for thin accretion discs 
where the vertical height limits the scale of the turbulence, the Prandlt law for viscosity works very well 



D Corrections to the radial temperature gradient 

Here we shall estimate how the radial temperature gradient differs from the adiabiatic law due to the 
convective motions in the shell. By multiplying ((56)) by (l/2)(a;m —a;*), we obtain the convective heating 
rate caused by interaction of the shell with the magnetosphere: 

L,^]^ZMR\{u:^-uj* f . (D.l) 

Multiplying the same (|56p with oj* yields the rate of change of the mechanical energy of the neutron star 

Lk^ ZMR\uj*{uj^-Lo*). (D.2) 

The total energy balance is then 

Lt^L, + Lu = ^ZMR\{loI - Lj*^) . (D.3) 

Note that the obtained formula for is similar to that describing energy release in the boundary layer 
of an accretion disc, see [79j. 
The convective energy flux is: 

_ ZMR\{u;^-^* f 

The convective energy flux can also be related to the entropy gradient as ^80] ) : 

qc = -^'^cT^ , (D.5) 
where S is the specific entropy (per gram). Here is the radial turbulent heat conductivity 

Vc =< uJc >= ChUcR , (D.6) 

where the characteristic scale of convection Ic ^ R, the velocity of convective motions Uc ~ Cs ^ R^^^^, 
and Ch is a numerical coefficient of the order of one. Thus 

= ^c(i?A) i^-j^j ■ (D.7) 
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Next, we make use of the thermodynamic identity for the specific enthalpy H: 



dH 



p dR 



dR ' 



(Di 



We remind the reader that the enthalpy can be written as 

dH = CpdT , 

where 

Cp 



T 



n 



9Tjp 7-1 Mm 

is the specific heat capacity at constant pressure. Expressing T{dS/dR) from (jD.5[) and making use of 
the hydrostatic equation [(U)] written as 



dPg/p 



n GM 



■0(7, TO*) • 



dR fJ-mCp R^ 

the thermodynamic identity (jP.Sp can be rewritten in the form 



dT 
dR 



1 



GM ZUr{RA)(RA 



i?2 



2iy,{RA) 



(D.9) 



By definition the adiabatic temperature gradient is determined by the first term on the right hand side 
{dT/dR)ad = g/cp. Equation ((DJt can be integrated to find the actual dependence of the temperature 
on the radius in the convective shell: 



GM , , , ZUr{RA) „3 



R 



2v,{R, 



Rj^{uJn 



R 

rI 



Close to equilibrium {lui* — 0), we can use (jB.12|) and write 



T = 



GM 



2GhU,{RA) 



*2 - 



In 



R 

rI 



(D.IO) 



(D.ll) 



This solution shows that in the whole region between Ra and Rb, for slowly rotating pulsars (i.e., in 
which LUjn <Si ljk{Ra))i the temperature distribution is close to the adiabatic law with a temperature 
gradient close to the adiabatic one [(0])]: 



T: 



7- 



IGM 



7 UR 



(D.12) 



Here we have only taken into account energy release due to the frequency difference near the magneto- 
sphere. In reality, there may be additional sources of energy in the shell (e.g. the heat release during 
magnetic reconnection and turbulence (see Appendix E), etc.). 



E Dynamics of a stationary spherically-symmetric ideal gas flow 

In this Appendix, we write down the gas-dynamic equations of a spherically symmetric ideal gas flow 
onto a Newtonian gravitating center. This problem was considered in the classical paper by Bondi f|18j) 
for an adibatic accretion flow. Adiabatic gas outflows (stellar winds) were studied by Parker [ST]. A 
thorough and comprehensible discussion of such flows can be found in the monograph by V. Beskin |82j . 
Here we focus on the role of the cooling/heating processes near the Alvenic surface, and also take into 
account the effects of turbulence and/or convection (anisotropy in general). As discussed in the main 
text, at low X-ray luminosities the quasi-static shell is capable of removing angular momentum from the 
rotating magnetosphere via convective motions. As the accretion rate exceeds some critical value, strong 
Compton cooling causes a free-fall gap to appear above the magnetosphere, and angular momentum 
cannot be transferred from the magnetosphere to the shell any more. 
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The equation of motion (|A.16p in the absence of viscosity reads: 

dur IdPa idPl 2(P*-Pl) GM 



dR 



p dR p dR 



pR 



i?2 



(E.l) 



Here Pg = pc^g/^ is the gas pressure, and P* the pressure due to turbulent pulsations, which in general 
are anisotropic: 

s'"|| (E-2) 



P|f = p < u\ pm\cl = lPgm\ 



Pi = 2p < u\ 2pmlc2 = 2jPgm^_ 



(E.3) 



< U( >=< m| > +2 < > is the turbulent velocity dispersion, m| and ml are the parallel and 
perpendicular turbulent Mach numbers squared). 
From the first law of thermodynamics we have 



dE _ Pg dp ^dS 
dR ^ ~p"dR ^ dR ' 



where the specific internal energy (per gram) is 

E = cvT 

and the heat capacity is 

cv 



7(7- 1) 

n 1 



(E.4) 



(E.5) 



(E.l 



Mm 7 - 1 

From the second law of thermodynamics, the specific entropy change can be written using the rate of the 
specific heat change dQ/dt [erg/s/g] as 



T 



dS _dQ _ dQ/dt 
dR dR Ur 



Using the mass continuity equation 
we find 



M = AirRyur , 



1 dp 
'pdR 



Using the relation = jT^T, we finally obtain: 



c2 dR 



(7-1) 



2 1 du^ 
R ^ 2ul dR 



2 1 dul 
R ^ 2^"Zr 



dQ/dt 

UrCyT 



(E.7) 

(E.8) 
(E.9) 

(E.IO) 



Note that this equation can also be derived directly from the ideal gas equation of state written in the 
form 



Pg = Ke^/"-' p''' , 

where K is some constant. 

Using (|ElO| . the gas pressure gradient can be rewritten in the form: 



1 dPg _ cl dQ/dt 

Pg dR CyUr T ^ 

Plugging (jE.12[) into the equation of motion finally yields: 



2 1 dul 
R ^ 2ul dR 



1 1 dul 
2^~dR 



c'(l + 7"^P 



2, / 2 dQ/dt\ 



R CyUrT J 



2cl 



{m\ - ml) GM 
R W 



/ ul- c^(l + 7mjj) 



(E.ll) 



(E.12) 



(E.13) 



Note also that in the strongly anisotropic case where m| — ml ^ ml , the role of turbulence increases 
in comparison with the isotropic case where m^ = ml = (l/3)ml 
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We can also introduce the Mach number in the flow M = Ur/cg. Then from (jE.10[) and (|E.13p we 
derive the equation for the Mach number: 



dR 

2[(7-l)A<''-(7+l)(mjj-mi)] _ [M^+yiM^^J^A _ (7+l)GM 
R cvT dR R^cl 



(E.14) 




where we have substituted {dQ/dt) — u{dQ/dR). Equations (jE.lOp . ([E.lSp and (IE. 141) can be used to 
solve the dynamics of the accretion flow for pairs of independent variables (m,Cs), {u,M.) or {cs,M). 
Here, however, we shall only consider the behaviour of the flux near the singular point. To this end, we 
can use (|E.13p . 

(IE.13[) has a singular saddle point where the denominator vanishes: 

ul = cl{l + ^m\). (E.15) 
So must the numerator, from which we find the quadratic equation for the velocity at the singular point: 

l + (7-l)™f+<\ (dQ/dt\ GM 

Remember that in the adiabatic case [dQ/dt — 0) without turbulence at the saddle point we have simply 

- c2 = — . (E.17) 

We stress that the presence of turbulence increases the velocity at the singular point. For example, for 
7 = 5/3 we find for strong anisotropic turbulence = Cg(l + (5/3)m|); for the isotropic turbulence the 
correction is smaller: v? — 0^(1 + (5/9)m(). The transition through the sound speed (the sonic point 
where = c^) lies above the saddle point due to turbulence, and there is no singularity in the sonic 
point. 

First let us determine the turbulence heating rate in the quasi-static shell {dQ/dt)f: 

dQV_l<^^,>^ (E.18) 



dt J ^ 2 tt 
where the characteristic time of the turbulent heating is 

R R 

tt^at— = at , (E.19) 

Ut mtCs 

with at being a dimensionless constant characterizing the turbulent dissipation energy rate and the 

= to| + 2r7i^ . The turbulent heating rate can thus be written as 



dQ 

dt If 2atR 



^= 3 

-nit 



In the case of Compton cooling we have 

' dQ\ cv{T — T,^) 



dt I ^ tc 



(E.20) 



(E.21) 




where tc is the Compton cooling time [ p7)) ]. 
(jE.16|) can now be written in the form: 

1 + (7 - l)m| + mj \ 7(7 - l)ml u,(l - T^/T) GM _ 

l + 7m| j ''ur 2atR -ftc R^ ' ^ ' ' 

As we study the accretion problem, the sign of the velocity Ur — dR/dt is negative, so below we shall 
write Ur = — |ur|- Then for the absolute value of the velocity at the singular point where the sound speed 
is Cs/|wr| = —1/(1 + 7"ipi/2 have the quadratic equation: 

ulU'^^'^ '^i^ "^^ V ' 2U/2 ' - '"-'^^ " "^'^"^^ ^ ^ - - (E.23) 

'^i? I l+7m2 i ^^(1 + 77712)1/2 2atR jtc R^ 
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In this case, the solution to (jE.16p reads: 



\Ur\ 



R{1 - T^/T) /2GM 



1 R R\l-T,/Tf 
AA 2GM m-fHj^A^ 



AjtcA V R 
where we have introduced the dimensionless factor 

1 + (7 - l)m| + m\ 7(7 - + 2m^)^/2 



1/2 



A = 



1 + 7TOy 



4at(l +7mj^)i/2 



(E.24) 



(E.25) 



In the case of isotropic turbulence where r7i|| = mj^ = 1/^/3, nit — 1, for 7 = 5/3 the factor A w 1.23, 
and in the case of strongly anisotropic turbulence where m|| = l,m± = 0,mt = 1, this factor is A « 0.8. 
In units of the free-fall velocity the solution (jE.24| reads: 



/(„) = ^ = (1 - T./T) /tjf^ 
Uff 4:jA \ tc 



1 ^ {l-TjTf (tff 



A 



tc 



1/2 



With Compton cooling present, the temperature changes exponentially: 



(E.26) 



(E.27) 



see the main text). When cooling is slow, t///tc ^ Ij the critical point lies inside the Alfven surface, i.e. 
no transition through the critical point occurs in the flow before it meets the magnetosphere, and in this 
case we expect settling accretion. If this point lies above the Alfven surface, the velocity of the flow may 
become supersonic above the magnetosphere, and one may thus expect the formation of a shock. Both 
turbulence and rapid cooling shifts the location of the critical point upwards in the flow. 

In the case of rapid cooling tfj/tc 1, T — >• T^,, so again Ur/uff « 1/2 (cf. (|E.17I) for an adiabatic 
flow), but the critical point now lies above the Alfven surface, so a free-fall gap above the magnetosphere 
appears. The ratio f{u) = \ur\/uff reaches a maximum at tfj/tc ~ 0.46 (assuming a typical ratio 
Tcr/Tx — 10), and depending on the value of A — 0.8-^ 1.23 (anisotropic or isotropic turbulence) it equals 
to /(u) = 0.5-0.6. 
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